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DESCRIPTION  OF  WORK:  The  objective  of  this  contract  is  to  develop  new 

robust  techniques  for  estimating  the  directions  of  arrival  of  multiple  signals 
utilizing  the  available  multi-sensor  information.  This  includes  direct  procedures 
utilizing  the  generalized  eigenvalues  associated  with  certain  matrices  generated 
from  the  signal  subspace  eigenvectors  of  the  actual  array  output  matrix.  In 
addition,  this  proposal  addresses  the  problem  of  analyzing  these  techniques  to 
evaluate  their  performance  when  the  array  output  cross-covariances  are  directly 
estimated  from  the  data.  In  this  regard,  the  mean  and  variance  expressions  for 
the  angle-of-arrival  estimators  can  be  used  to  derive  thresholds  for  resolving  two 
or  more  closely  spaced  sources. 

Further,  by  interpreting  the  maximum  entropy  method  (MEM)  of  spectrum 
estimation  in  terms  of  maximization  of  minimum  mean  square  error  for  a  one 
step  predictor,  a  new  class  of  spectrum  extension  problems  to  identify  the  best  r- 
step  predictors  is  proposed. 

FY88  PROGRESS:  A  comprehensive  asymptotic  analysis  of  a  class  of  high 

resolution  estimators  for  resolving  correlated  and  coherent  sources  in  white  noise 
has  been  completed.  Using  these  results  resolution  thresholds  have  been  derived 
in  two  and  three  source  scenes  for  uncorrelated  as  well  as  coherent  situations. 

Given  a  finite  set  of  n  autocorrelations  of  a  stationary  discrete-time 
stochastic  process,  the  well  known  problem  of  extending  this  given  sequence  so 
that  the  power  spectral  density  associated  with  the  resulting  infinite  sequence  is 
nonnegative  everywhere  is  further  investigated.  Motivated  by  the  maximum 
entropy  extension,  which  is  equivalent  to  the  maximization  of  the  minimum 
mean  square  error  (MMSE)  associated  with  one-step  predictors,  the  natural 
extension  of  maximizing  the  MMSE  for  r-step  predictors  compatible  with  the 
given  correlations  is  studied.  This  analysis  shows  the  resulting  spectrum 
corresponds  to  that  of  a  stable  ARMA  (n,  r  -  1)  process. 
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Technical  Report 


This  research  report  concerns  primarily  with  new  algorithms  in 
estimating  signals  from  sensor  output  measurements.  A  complete  performance 
analysis  of  the  proposed  algorithms  is  developed  in  Appendix  A  together  with  new 
results  for  resolving  two  and  three  uncorrelated/coherent  signals. 

Appendix  B  concerns  with  the  problem  of  extending  a  given  set  of  n 
autocorrelations  so  that  the  power  spectral  density  associated  with  the  resulting 
infinite  sequence  is  nonnegative  everywhere.  Motivated  by  the  maximum  entropy 
method,  which  is  equivalent  to  the  maximization  of  the  minimum  mean  square 
error  (MMSE)  associated  with  one-step  predictors,  the  natural  extension  of 
maximizing  the  MMSE  for  r-step  predictors  compatible  with  the  given 
correlations  is  studied.  It  is  shown  here  for  the  first  time  that  the  resulting 
spectrum  corresponds  to  that  of  a  stable  ARMA  (n,  r  -  1)  process.  ^The  details  of 
this  optimum  filter  for  a  two  step  predictor  are  presented  in  Appendix  B  along 


with  several  other  interesting  conclusions. 
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Appendix  A 

New  Resolution  Threshold  Results  in 
Three  Source  Scenes 


Abstract 

This  paper  presents  new  performance  analysis  results  in  a  three-source  scene,  when 
MUSIC-type  high  resolution  estimators  are  used,  to  estimate  the  directions  of  arrival  of 
incoming  signals.  A  two-source  scenario  has  been  popularly  used  to  measure  the  quality 
of  various  estimators  in  terms  of  the  SNR  required  to  resolve  {resolution  threshold)  two 
closely  spaced  sources.  Similar  results  in  a  more  realistic  three-source  scene  is  presented 
here  for  uncorrelated  as  well  as  coherent  situations  along  with  some  interesting  possibili¬ 
ties  that  deserve  further  study. 


I.  Introduction 


Given  a  set  of  alternatives  that  perform  identically  under  ideal  conditions,  one 
important  question  is  how  to  measure  their  robustness  or  sensitiveness  to  imperfec¬ 
tions  that  will  invariably  exist  in  reality.  To  illustrate  this,  consider  the  array  process¬ 
ing  scene  where  a  set  of  sensors  collect  data  from  signals  present  in  their  field  of  view 
in  presence  of  noise.  Under  the  assumption  that  the  sensor  noise  are  independent  and 
identical  between  themselves,  a  variety  of  high  resolution  techniques  have  been 
developed  in  the  recent  past  [1  -  5]  to  estimate  the  directions  of  arrival  of  incoming 
signals.  They  all  depend  on  the  fact  that  the  covariance  matrix  formed  from  the  sen¬ 
sor  array  output  data  has  several  interesting  structural  properties.  For  example,  in  an 
equal  noise  situation  when  none  of  the  signals  present  in  the  scene  are  completely 
coherent  with  each  other,  the  direction  vectors  associated  with  the  actual  arrival 
angles  can  be  shown  to  be  orthogonal  to  the  eigenvectors  corresponding  to  the  lowest 
eigenvalue  of  the  array  output  covariance  matrix.  This  forms  the  basis  for  the  MUSIC 

*  The  authors  are  with  the  department  of  Electrical  Engineering  and  Computer  Science  at 
Polytechnic  University,  Brooklyn,  New  York.  This  work  was  supported  by  the  Office  of  Naval 
Research  under  contract  N-00014-89-J-1512. 
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algorithm  [1]  and  several  of  its  successors. 

The  coherent  source  situation  can  also  be  brought  under  this  formulation,  by 
employing  certain  preprocessing  on  sections  of  the  array  output  covariance  matrix.  In 
this  regard,  the  subarray  averaging  scheme  [6-8]  effectively  creates  a  smoothed  covari¬ 
ance  matrix  that  is  structurally  identical  to  a  covariance  matrix  in  some  correlated 
situation.  Even  more  resent  schemes  such  as  ESPRIT,  TLS-ESPRIT  [3,  9],  that  make 
use  of  an  underlying  rotational  invariance  property  of  certain  vectors  in  the  signal  sub¬ 
space  stems  from  the  above  idea.  More  interestingly,  under  ideal  conditions,  all  these 
techniques  perform  identically,  i.e.,  when  the  ensemble  averages  of  the  array  output 
covariances  are  exactly  known,  any  source  situation  that  is  resolvable  by  any  one  of 
these  techniques  can  also  be  resolved  by  any  other  scheme  that  belongs  to  this  general 
category,  in  contrast  to  -  say  -  the  linear  prediction  method. 

However,  the  situation  is  entirely  different  when  array  output  date  is  used  to  esti¬ 
mate  these  covariances.  In  that  case,  all  these  techniques  can  be  viewed  as  distinct 
algorithms  and  they  will  perform  differently  for  the  same  source  scene.  Several  cri¬ 
teria  have  been  proposed  to  evaluate  their  performance  under  such  conditions.  Of 
these  the  bias  and  variance  of  the  estimator  give  reasonable  information  about  the 
robustness  of  the  estimator  under  consideration  compared  to  other  techniques.  A 
more  meaningful  physical  measure  under  these  circumstances  is  the  signal-to-noise 
ratio  (SNR)  required  to  resolve  two  closely  spaced  sources.  At  least  qualitatively,  this 
information  tells  about  the  sensitivity  or  superiority  of  the  technique  under  considera¬ 
tion.  The  bias  and  variance  expressions  described  above  has  been  used  in  this  connec¬ 
tion  to  evaluate  the  resolution  threshold  in  a  variety  of  two  source  scenarios  [  10-13]. 

Naturally,  two  source  scene  is  often  unrealistic,  and  to  evaluate  the  degradation 
in  performance  of  these  estimators  in  a  general  set  up  it  is  necessary  to  analyze  at 
least  a  three-source  scene.  Reevaluating  the  resolution  threshold  in  such  a  scene  will 
also  indicate  the  relative  degradation  in  performance  because  of  the  presence  of  yet 
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another  source  and  in  this  paper  we  propose  to  address  this  problem  in  uncorreiated 
and  coherent  situations. 

II.  Main  Results  in  a  Three-Source  Scene 

In  this  section  we  will  consider  three  source  scenarios  in  two  extreme  situations, 
viz,  the  uncorreiated  case  and  the  coherent  case.  The  uncorrelated  case  represents  the 
best  possible  scene  in  terms  of  the  SNR  required  to  resolve  three  signals  and  the 
coherent  scene  represents  the  worst  scenario  because  of  the  complete  dependence 
among  the  signals.  As  remarked  earlier,  the  coherent  scene  can  be  decorrelated  by 
employing  the  forward/backward  smoothing  scheme  on  the  subarray  output  covari¬ 
ance  matrices  [8].  Though,  in  general,  this  requires  additional  number  of  sensor  ele¬ 
ments,  interestingly  any  source  scene  containing  two  coherent  signals  can  be  decorre¬ 
lated  without  the  use  of  extra  sensor  elements.  In  fact,  the  mean  of  the  array  output 
covariance  matrices  corresponding  to  the  forward  array  and  its  complex  conjugated 
backward  version  can  always  resolve  any  such  source  scene  [8]  and  for  uniformity  in 
comparison,  we  will  adopt  this  situation  in  this  study.  Thus,  the  coherent  scene  under 
discussion  here  will  consist  of  two  coherent  signals  that  are  uncorreiated  with  the  third 
signal.  The  analysis  that  follows  in  that  case  will  utilize  the  forward/backward 
smoothing  scheme  that  employs  the  forward  array  and  its  complex  conjugated  back¬ 
ward  version  to  decorrelate  and  resolve  the  three  signals  and  their  arrival  angles.  To 
begin  with,  we  consider  the  uncorreiated  source  scene. 

IIA  Uncorreiated  Signals 

Consider  three  uncorreiated,  planar  source  waveforms  ),  u  2(t )  and  u3(t ) 
arriving  at  an  M -element  uniform  linear  array  along  the  directions  $v  $2  and  03  with 
respect  to  the  line  of  the  array.  Under  i.i.d.  noise  conditions  the  M  xM  array  output 
covariance  matrix  has  the  form(1)  [8] 


(1) 


R  =  ARwAf  +  ol 

whe„-  the  uncorrelated  condition  among  the  sources  implies  that,  the  source 
covariance  matrix  has  the  form  Rw  =  diag  [Pv  P2,  P3  ]  where  P  =  E  [  f  uj  (t )  |  2],  i  = 
1,  2, 3  represents  the  signal  power.  Moreover 

A  =  VM  [  a^),  a(u/2),  a(w3)  ]  (2) 

where 


a  (w*)  = 


Vm  ' 


[l,e 


•jvk  -j  2uk 


-;( Af-l)w*  T 

e  ] 


UJ,. 


-  7T  COS  0, 


(3) 


represents  the  direction  vector  associated  with  the  arrival  angle  Qk  and  o~  = 
E  [  |  n:  (t )  |  ],i  =  1, 2,  •  •  • ,  M  denotes  the  common  noise  variance. 

If  >  ^  >  •  •  •  >  \M  and  0V  02,  ,  0M  represent  the  eigenvalues  and 

corresponding  normalized  eigenvectors  of  R,  i.e., 

R  -  s  A.  fi,  f>)  (4) 

/-l 

2  t 

then  in  this  particular  setup  \  -  a  for  i  >3  and  moreover  0\  a (uk )  =  0,  for  all  i  >  3 
and  k  =  1,  2, 3.  Thus,  the  zeros  of 


M 


M  *  *>  J  f 

CM-  E  1  fij  aM  1 2  =  i  ~  E  l^aM 


(5) 


represents  the  true  angles-of-arrival  [1]  and  in  principle,  this  estimator  can  resolve  any 
three  sources  irrespective  of  their  angular  separation. 

However,  when  array  output  data  sample  vectors  are  used  in  estimating  R  by 
means  of  any  standard  procedure,  the  corresponding  estimator  in  (5)  is  a  random  vari¬ 
able  and  the  above  conclusions  are  only  approximately  true.  In  fact,  when  the 

(1)  Hereonwards  a  (or  A),  a  and  A  stand  for  scalar,  vector  and  matrix  in  that  order. 
Similarly,  A*,  A  and  Af  represent  the  complex  conjugate,  transpose  and  complex  conjugate 
transpose  of  A  respectively. 
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maximum  likelihood  method  is  used  to  estimate  the  covariance  matrix  R  using  N 
data  vectors,  the  associated  estimates  of  eigenvectors  so  obtained  can  be  used  in  (5)  to 
generate  the  sample  estimator  The  statistical  properties  of  Q(u>)  has  been  well 
documented  in  the  case  of  zero  mean  complex  circular  Gaussian  data  for  the  general 
setup  when  smoothing  is  employed  on  the  subarray  covariance  matrices  to  decorrelate 
coherent  signals  present  in  the  data  [11,  12].  As  a  special  case  of  this  general  analysis, 
the  bias  and  variance  of  the  sample  MUSIC  estimator,  Q  (w)  has  been  shown  to  be 
[11] 

A  <72 

£[(2(u)]  =  Q(u/)  +  4s  - {(M  -K)l0>a(u)l2-Q(<j)}+O(-JT)  (6) 

;-l  (X  -af  ' 

and 

A  a 

Varf&m  =  £ - L—  |j9/a(u,)|2  +  0(-4).  (7) 

Thus,  in  particular,  within  the  above  first-order  approximation  Var(Q(uk))  =  0,  1  < 
k<  K.  However,  along  the  true  angles  of  arrival,  E[Q{uk)\  f  0  (see  (6))  and  the 
deviation  of  E  [Q  (uk )]  from  zero  -  their  nominal  value  -  suggests  the  loss  in  resolution 
for  this  estimator  below  a  certain  angular  separation.  Since  the  estimator  has  zero 
variance  along  the  true  arrival  angles,  for  a  fixed  number  of  samples  a  threshold  in 
terms  of  SNR  exists  below  which  the  nulls  corresponding  to  the  true  arrival  angles  are 
no  longer  separately  identifiable.  Considering  two  sources  at  a  time,  the  correspond¬ 
ing  sources  are  separately  identifiable  if  the  bias  at  their  middle  angle  is  larger  than 
the  maximum  of  that  at  either  of  the  two  arrival  angles.  Letting  9X  <  02  <  0y  this  gives 
rise  to  the  following  inequalities  for  resolving  three  sources: 

£[^((^1  +  w2)/2  )  ]  >  max  {E  [  @  (u/j)  ],E[Q  (o/2)  ]}  (8) 

and 
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(9) 


E  [  £  ( (uj  +  uj)/2 )  ]  >  max  {£  [  (w2)  ] ,  £  [  Q («,)  ]>  ■ 

When  the  sources  are  equispaced,  i.e.,  w1  -  u/2  =  -  w3,  the  minimum  value  of  SNR 

that  satisfies  both  (8)  and  (9)  may  be  taken  as  the  resolution  threshold  in  a  three- 
source  scene.  In  the  special  case,  when  all  signal  powers  are  equal  (  Pi  =  P,i  =1,  2, 
3),  from  the  resulting  symmetry  it  can  be  established  that 

£[2(H  +  «2>/2)]  =E[Q((u2  +  «j)/2)], 

and  (8)  and  (9)  collapses  into 

£t!2((«t+“2)/2)]>£[^("2))-  <10> 

2 

Hence,  by  definition,  the  desired  (normalized)  resolution  threshold  £  =  MP  jo 
satisfies  (10)  with  an  equality  and  in  that  case  one  can  resolve  three  equipowered, 
uncorrelated  sources. 

To  complete  the  analysis,  from  (6),  it  remains  to  obtain  explicit  expressions  for 
the  eigenvalues  and  eigenvectors  associated  with  the  signal  subspace  of  R  and  we 
proceed  to  do  so  in  the  next  section. 

II.B  Eigenparameters  in  a  Three  Uncorrelated  Source  Scene 

To  begin  with,  notice  that  (from  (1))  the  three  largest  eigenvalues  A^  A,  and  A.  of 
R  are  related  to  the  eigenvalues  ^  and  ^  of  the  3x3  matrix,  Ru  ArA  through  the 
relation 

Af  =  fij  +o2  ,  i  =  1,2,3.  (11) 

In  an  equipowered  situation,  however, 


Ru  A1  A  =  MP 


1  P 12  ^13 

r 

0  ^12  ^13 

\ 

Pl2  1  p73 

-  MP 

I3  + 

^12  0  p2 3 

P\3  P23  1 

V 

P\3  p22  0 

(12) 


where 


Pik  -at(w|.)a(wJk)  ,  i,k  =1,2,3 

represents  the  spatial  correlation  coefficient  between  the  Ith  and  k,h  sources.  Thus,  if 
uv  u2  and  represent  the  eigenvalues  of  the  zero  axial  matrix 


0  ^12  ^13 
P\2  0  Pn  ’ 


P13  P73  0 


(13) 


then  from  (12) 


Hi  =  MP(  1  +  v.)  ,  /=  1,2,3. 


(14) 


Towards  obtaining  the  eigenvalues  of  the  above  zero  axial  matrix,  observe  that  in  an 
equispaced  source  situation  (i.e.,  -  u2  =  u>2  -  u3  -  2ud ), 


p12  ~  P73  ~  e  Ps 


sinA/cj. 

,  4 - i 

5  M  sinu^ 


(15) 


and 


j  2(M-l)wd  sin2^Wd  A  COsMud 


COS  U, 


(16) 


With  (15)  and  (16)  in  (13),  its  eigenvalues  satisfy  the  cubic  equation 
|  B  -  i/I3  j  =  v  -  p]{2  +  p])v  -  2 pispc  =  0. 
Interestingly,  the  above  equation  can  be  factored  into  the  form 

iy  +  ps pc) (y  -  ps pc u~  2Ps) 

and  hence  the  three  real  roots  of  (17)  are  given  by 


(17) 


"2  =  -Ps  Pc 


-11- 


which  gives  the  desired  eigenvalues  in  (14)  to  be  (in  decreasing  magnitude  when  ps ,  p: 
are  positive;  i.e.,  for  small  angular  separations) 


h2=MP(1  —ps  pc ) 


(lS.a) 

(18.b) 

(18.c) 


From  (4)  and  the  discussion  that  follows,  it  is  clear  that  the  eigenvectors  associated 
with  (11)  are  linear  combinations  of  the  true  direction  vectors,  Le., 

Pi  «  aK)  +  ku  a(w2)  +  k^  a(w3)  ,  i  =  1,  2,  3  (19) 

and 

(AR«  A')0i  m  M  p(a(wi)  at(wi)  +  a(w2)  aV2)  +  a(w3)  af(o;3) )  0.  =  ^  /?■  i  =  1,  2,  3 . 

(20) 

With  (19)  in  (20)  and  equating  coefficients  of  the  true  direction  vectors,  we  obtain 
three  consistent  equations  for  each  i .  Using  two  of  these  equations  the  unknowns  k  u 
and  k^  in  (19)  can  be  evaluated  (for  details,  see  section  II.C)  and  after  some  algebra, 
this  gives 

.j(M- l)Wrf  ) 

P2  «  [e  a(^i)  -  e  a(u/3)J 


-12- 


Letting 


j(M-  Du%/2  ...  , 

u;  =  e  a(w- )  i  =  1, 2, 3 


(21) 


and  normalizing  the  above  vectors,  with  the  help  of  =  uju3  =  ps,  u[u3  =  ps  pc, 
we  obtain  the  desired  eigenvectors  to  be 


Ul"U3 


A- 


y/2(l-PsPc) 


(22) 


and 


Cl/  (U1  +  u3>  +  C  2i  **2 
y/  '2 Cl  (!  +  Ps  Pc)  +  40 li; C2;  Ps  +  C2 


t  =1,3 


where 


(23) 


Pc2  PcV^i  3pc  >/  Pc  +  8 

cii  =  1  +  T  * - 2 -  5  c2/  =  T"  *  - 2 - • 

Notice  that  /?2  is  orthogonal  to  u2  and  is  a  linear  combination  of  u:  and  u3  only. 

As  remarked  earlier,  the  desired  resolution  threshold  satisfies  ( 10)  with  an  equal¬ 
ity  and  in  an  equipowered,  equispaced  three  source  scene  the  above  eigenparameters 
can  be  used  to  evaluate  that  explicitly.  Towards  this  purpose  using  (11)  and  (14)  we 
get 


,\  2\2 


(/*,-  +  o)° 

2 

Hi 


1 


e2/,2 


(24.a) 


where  the  array  output  signal-to-noise  ratio  is 

£  ±MP/a 


and 


/,•  -  1  +  v.  ,  /  =  1, 2, 3 . 


(24. b) 
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Thus  with  A-  =  3,  from  (6) 


E&M]  = 


M±1  f 

N  S 


(  \p]&(u2)  | 2  |  a(wj)  | 


f',- 


i2!,1 


*0lW) 


(25) 


and  similarly  with  um  =  (w1  +  w2)/2  representing  the  middle  angle,  we  have 


£[(?(«„)] -CK)  +  iri: 


f  (M-3)|^a(Um)|2-Q(K„)) 


1=1 


(M  -3)  |  pj a(um  )\  -Q  (ujm ) 

?? 


+  0(^ 


(26) 


Finally,  equating  (25)  and  (26)  and  neglecting  terms  of  order  less  than  1  /N ,  we  obtain 
the  quadratic  equation 

a,t f  +  63f  +  c3  =  0  (27) 

and  this  gives  the  desired  resolution  threshold  output  SNR  in  an  equipowered, 
equispaced,  uncorrelated  source  scene  to  be 

-  b*+  Vb}  -  4a, c, 

«3  = - 25; - •  <28> 


Here 

.  I  3  (jW-3)(|^/a(um)|2-  |  0*  a(«j)  I  2)  ) 

*34jvS - r - 

1*1  1 

\ 

A  1  ±  (Af  -3)(  |  )  | 2  -  I  ^/a(w2)  | 2)  -Q(ullt ) 

C3  a  W  E  P 

1=1  li 
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and  /,. ,  0. ,  i  =  1, 2, 3  are  as  in  (22)-(24). 

The  corresponding  threshold  results  in  a  two  source  scene  can  be  utilized  to  esti¬ 
mate  the  degradation  in  performance  created  by  the  presence  of  a  third  source.  In  the 
case  of  two  equipowered  uncorrelated  sources,  the  resolution  threshold  has  been 
shown  to  be  [10], 

1  f  20(Af  -  2) 

N\  a4 

where  •* 

A2  =  A/2w2/3  . 

Though  and  f2  possess  similar  features,  for  a  given  array  the  resolution  threshold  in 
the  three  source  case  can  be  much  larger  than  that  in  the  two  source  case.  The  above 
asymptotic  analysis  is  also  found  to  be  in  agreement  with  Monte  Carlo  simulation 
results.  Table  1  represents  a  typical  case  study.  Simulation  results  indicate  that  when 
equality  holds  in  (10),  the  probability  of  resolution  ranges  between  33  to  50  percent  in 
both  cases.  This  suggests  that  the  above  results  should  give  an  approximate  threshold 
in  terms  of  £  for  0.33  to  0.5  probability  of  resolution  region  and  comparisons  are  car¬ 
ried  out  in  Fig.l  using  (28),  (29)  together  with  simulation  results  from  Table  1  for  this 
range  of  probability  of  resolution.  Fig.  2  shows  another  comparison  for  a  different 
array  length.  From  these  data  it  is  also  clear  that,  for  the  same  SNR,  sources  that  are 
located  around  the  broad  side  of  the  array  can  be  much  closer  than  those  arriving 
along  the  line  of  the  array.  Notice  that  in  all  these  comparisons,  the  source  locations 
have  been  chosen  such  that  they  are  always  well  within'the  main  beam  of  the  array. 

Interestingly,  these  results  indicate  that  the  presence  of  a  third  source,  that  is 
symmetrically  located  with  respect  to  the  center  source,  increases  the  SNR  required  to 
resolve  the  original  two  sources  approximately  by  a  factor  of  two  (in  dBs).  This  in 
turn  speaks  for  the  sensitivity  of  the  MUSIC-type  eigenstructure  based  algorithms. 


1  + 


N 


5(M  -  2) 


(29) 
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In  the  next  subsection  we  proceed  to  analyze  the  coherent  case,  where  two  of  the 
three  signals  present  in  the  scene  are  completely  coherent  with  each  other  and  the 
third  signal  is  uncorrelated  with  these  two  signals. 

II.C  Two  Coherent  Signals  (and  an  Uncorrelated  Signal) 

As  remarked  earlier,  the  coherent  scene  under  discussion  here  consists  of  two 
coherent  signals  that  are  uncorrelated  with  the  third  signal.  The  coherent  signals  can 
be  the  outermost  ones  or  the  adjacent  ones,  and  since  only  one  of  these  situations  is 
symmetric,  these  two  cases  have  to  be  analyzed  separately.  Thus  (a)  u^t)  and  u3(t) 
are  coherent  and  u2(t )  is  uncorrelated  (the  symmetric  case)  or  (b)  u  3(t )  and  u2(t )  are 
coherent  and  u3(t)  is  uncorrelated  with  the  other  two  signals.  Of  these  two  situations, 
we  will  only  deal  with  the  seemingly  ’simpler’  symmetric  case,  the  details  of  which 
itself  turn  out  to  be  quite  involved. 

To  begin  with,  in  the  symmetric  case,  consider  the  perfectly  coherent  situation, 
i.e.,  u3(t)  =  u x(/ )  and  the  mid-signal  u?(t)  is  uncorrelated  with  the  adjacent  ones.  By 
employing  the  covariance  matrices  corresponding  to  the  forward  array  and  its  complex 
conjugated  backward  version,  the  above  coherent  situation  can  be  decorrelated.  In 
that  case  the  smoothed  covariance  matrices  associated  with  the  mean  of  the  above 
covariance  matrix  has  the  form  [12] 


R  = 


R;  +  Rb 


=  P  A 


1  Op, 

0  1  0 


•  * 


P,  0  1 


Af  +  a\ 


(30) 


where 


-j (M  - 1) W3 

1  +  e  e 


2 


(31) 
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represents  the  effective  correlation  coefficient  due  to  f/b  smoothing  Here  R^  and 
R6  stand  for  the  forward  and  backward  array  output  covariance  matrices.  Clearly, 
is  the  same  as  R  in  (1)  and  with  $  as  in  (A.9) 

Using  (2)  and  (31)  in  (30),  we  have 

R=  MP  (a^) at(w1)  +  p' a(w3) a^i)  +  P,  a(wx) at(o,-3)  +  a(w3) af(u>3)  +  a(w2) af(w2) ]  +  a1 1 
=  b3bj  +  t^b*  +  b3b3)  +  a2I  =  BBf  +  <j2I 

where 

b1  -  a(wj)  +  a(w3) 
b2  ^v7a(w2) 

A  i(¥-  1)W1  /  x  .  /(M-l)w, 

b3  =e  a(wx)  +  e  a(w3) 

and 

M  .  _  ,t  .  2 

Since  R  =  2  \  with  \  -  o  for  i  >  3,  again,  the  signal  subspace  eigen- 

i-i 

values  Alt  A^  A3  of  R  are  related  to  the  three  positive  eigenvalues  jiv  jX,  £3  of  Bf  B/2 
through  the  relation, 

\  •MPji.  +  cr2  ,  i  =  1,2,3.  (32) 

However, 

(2)  Clearly,  | p,  |  should  not  be  equal  to  1,  or  the  angular  separation(=  2ud)  should  not  be  in 
the  neighborhood  of  kir/(M-l),  k  *  0,  1,  2,  ••• .  However,  in  general  u3(t)  =  a  ux(t), 
where  a  represents  the  complex  attenuation  and  in  that  case  the  above  restriction  does  not 
hold. 
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PU  P12  P13 


i  u  - 


®  ®  =  p12  p22  P23 

Pl3  p23  hi 


where  by  definition 


p..  =  ^-b^b. 
nj  2  t  j  ’ 


and  in  particular  with 


P,  -  cos  2  (M  - 1)  ujd  , 


and  ps,pc  as  defined  before,  it  is  easy  to  show  that 

hi  *  P33  *  1  +  ps  pc  pt 


(34.a) 


hi  ~  ^Ps  cos  W  ~  l)^d 


.  _  j(M- iyut_ 

P23  ~  e  P\2 


P\3  =e  ( PsPc+Pt )• 


It  now  follows  that  the  eigenvalues  of  B*  B/2  can  be  written  as 

h  =  ^11  -  °i  =  (1  +  Ps  Pc  Pt)  ~  °i  .  *  =  1. 2, 3 
where  */•  are  given  by  the  real  roots  of  the  cubic  equation 


(34.b) 


(34.c) 


(34.d) 


(34.e) 


U  p\2  P\3 

P\2  °~PsPcPt  p23  = 

P\3  P'23  V 
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V  -  PsPcPt(y  -  I  Pl3  I  2)  ■  (2  I  P]2  I  2  +  I  Pl3  I  V  +  2Re(PnPl3^  =  °*  (36) 
With  (34.a)  -  (34.e)  in  (36),  it  reduces  to 

V  -  PsPcPti*2  -  ( PsPc  +  Pt)2]  -  (2^2(X  +  Pt)  +  (PsPc  +  Pt)2)  0  +  2Ps(l  +  P')(PsPc  +  Pi) 

■  (*  ”  (P5PC  +  P,)}  +  ((p5Pc  +  Pf)  “  PsPcPt]  0  ~  2Ps(l  +  Pt)  ~  PsPcPMPc  +  P,))  =  °- 

Thus,  the  roots  of  the  above  equation  are  given  by 

^  °1  ~  Ps  Pc  +  Pt  (37-a) 

*1,3  =  "  T  ( <*/»c  +  Pt  ~  PsPcPt)  ±  '/D  )  (37-b) 

where 

D  =  (p*Pc  +  Pt  +  PsPcPt  )2  +  8ps2(  1  +  Pt )  •  (37.c) 

Finally  with  (37.a)  -  (37.b)  in  (35)  and  (32),  we  get  the  desired  signal  subspace  eigen¬ 
values  to  be 

Aj  -MP  Aj  +  a1  -  MP  (1  -  ^>(1  -  />,)  +  <r2  (38) 

Alt3  =  M  P  +  o’  ~  ~~~  (  1  +  (1  +  pspc)(\  +  pt)  ±  '/D  j  +  T2  (39) 

with  D  as  defined  in  (37.c). 

Towards  obtaining  the  eigenvectors  0V  02,  associated  with  these  eigenvalues, 
notice  that  once  again  they  are  linear  combinations  of  the  true  direction  vectors  or 
equivalently  those  of  uv  u^,  U3  with  u;  as  given  by  (21).  Thus, 

Pi  oc  ux  +  k^j  u3  ,  /  =  1,  2,  3 .  (40) 

To  make  further  progress,  we  rewrite  the  smoothed  covariance  matrix  R  also  in  terms 
of  ut,  iu  and  U3.  Using  (21)  and  (33)  in  (30),  this  gives 
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R  =  M ?[  (Uj  +  p,  u3)uj  +  +  (p,  ux  +  u3)u3t]  +  a\.  (41) 

Since 

R h  =A jh  ,  i  =  1,2,3,  (42) 

with  (32)  and  (41)  in  (42),  it  reduces  to 

( (uj  +  p,  H3>uJ  +  v^ul  +  ( p,  ux  +  u^)  /?.  =  /xf  0i  ,  i  =  1,  2,  3. 

Finally,  using  (40)  and  equating  coefficients  of  u.  on  both  sides  of  the  above  equation, 

Jr 

we  obtain 


P,(1+P,)  PSPc+P, 

k2 i‘ 

£/-( 1  +  PsPcPt) 

l~Pi  Ps 

= 

~Ps 

Ps(1+P t)  1+PsPcPt~Pi 

■k2i  ■ 

-( PSPc  +  Pt) 

For  reasons  that  will  become  apparent  soon,  using  the  first  and  third  row  this  gives 


k2 i 

((P,PC  +  Pf)2  -  O  +  PsPcP,  -  A-)2)/A/ 

-V 

1 

(44) 


provided 

A(.  =  ps  (1  +  p, )  ((1  -  pspc)(  1  -  pf)  -  fit)  7*0. 

Since  Aj  =  0  and  Ax  and  Aj  are  nonzero,  using  (40)  and  (44),  the  eigenvectors  0y  and 
03  can  be  expressed  as 

0i  <xps  (1  +  pt)  [(1  -  p5pc)(1  -  pt )  -  ) )  («!  +  u3) 
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-  ((!  -  PsPc)(l  ~  Pt)  ~  A-)  C1  +  PsPcPt  -  A  +  PsPc  +  Pi  )U2 

<xps  (1  +  pt){ 14  +  U3)  +  [ fi;  -  (1  +  Pspc)(  1  +  P,))u2  for  i  =  1,  3. 

By  making  use  of  the  exact  expressions  for  3  from  (39)  in  the  above  expres¬ 
sion,  the  corresponding  normalized  eigenvectors  take  the  form 

Ps  (1  +  P,)(u  1  +  u3)  +  A  “2 

A  = -  (45) 

11 A  11 

where 

1/2 

Cj-  =  J  (  1  -  (1  +  P^c)(l  +  P,)  ±  (  (1  -  (1  +  PsPc'X1  +  Pf)  )2  +  8^2(1  +  P;))  ) 

(46.a) 

and 

||  A  II  -  2p/(l  +  p,)2(l  +  PjPc)  +  4p2 (1  +  pf)di  +  c-2.  (46.b) 

Next  we  turn  our  attention  to  evaluate  ^2.  Since  the  determinant  corresponding 
to  the  first  two  equations  in  (43)  is  not  zero,  they  can  be  utilized  for  this  purpose.  In 
that  case,  proceeding  as  before  it  is  easy  to  show  that 

^2ocfc1u1  +  k2n2  +  ^3113 

where 

kx  =  pf(  1  +  pt)  +  (i \  -  l)(pspc  +  Pt ) 

k3  =  (P 2  -  l)2  -  PsPcpM 2  "  !)  “  P,2(!  +  Pt )  • 

By  making  use  of  (38),  it  follows  that  =  0  and  k3  =  -  fe  which  gives 

02aui  _  U3 
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or,  finally 


U1 

V2(l-Pt0c) 


(47) 


Clearly  (38)  -  (39)  together  with  (45)  -  (47)  completely  characterizes  the  signal  sub¬ 
space  eigenparameters  in  a  three-source  scene  consisting  of  two  symmetrically  located 
coherent  signals  and  a  centrally  located  uncorrelated  signal  in  terms  of  the  input  signal 
specifications  and  the  array  geometry. 

The  bias  and  variance  expressions  in  the  f/b  case  can  be  used  to  obtain  the 
corresponding  resolution  threshold  in  this  case.  General  expression  for  bias  in  the  f/b 
case  has  been  shown  to  be  [11,12] 
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and 
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where 


<2M-l-£  l#a(u)|2, 
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(50) 


and  is  the  inverted  0*  vector  with  *j.  m  -  0*M _m  +v  As  shown  in  Appendix  A,  7. , 
i  =  1,  2,  •  •  • ,  M  form  a  complete  set  of  eigenvectors  for  R  and  using  this  together 
with  certain  relations  regarding  equivalence  of  eigenvector  sets  of  a  matrix,  for  the 
above  source  scenario  it  is  shown  in  Appendix  B  that 
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where 

4P2Af2sin22(M  -  l)wd 

5(w)  =  -  |  ^2ui  I  2  Re  [ J(u)0zfa{u)]  . 

(X.-A,)^-^) 

2 

Notice  that  the  additional  term  B(u)  in  (51)  is  independent  of  the  noise  variance  a". 
Once  again,  following  (24.a)  we  obtain 
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Finally,  equating  (53)  and  (54)  and  rearranging  the  terms,  the  resolution  threshold 
associated  with  a  three-source  scene  consisting  of  two  symmetrically  located  coherent 
signals  and  an  uncorrelated  center  signal  satisfies  the  relation 

+  63e  +  c3  =  0.  (55) 


Here 
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The  signal  subspace  eigenparameters  that  have  been  obtained  in  (38)  -  (39)  and  (45) 
-  (47)  can  be  used  to  compute  the  desired  resolution  threshold  in  this  case  exactly. 
Once  again  the  corresponding  threshold  expression  in  a  two-coherent  source  scene 
can  be  utilized  to  estimate  the  degradation  in  performance  by  the  presence  of  an 
uncorrelated,  third  center  signal.  In  the  case  of  two  coherent  sources,  the  resolution 
threshold  has  been  shown  to  be  [11] 


N 


20  (Af  -  2) 

r  i 

_L_  A2' 

1  + 

’  N{M  -  2)  ' 

\ 

A4 

*  3A2 

20  16. 

5 M{M  -  4)  , 

/ 

.  (57) 


Table  2  represents  a  typical  case  study  for  a  three-source  symmetric  coherent 


scene  described  above.  From  these  simulation  results,  equality  between  (53)  and  (54) 
corresponds  to  0.33  to  0.5  probability  of  resolution  and  Fig.  3  shows  comparisons 
between  (55)  and  (57)  for  this  probability  of  resolution. 
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Equation  (55)  is  derived  based  on  the  assumption  that  inequalities  in  (8)  -  (10) 
are  satisfied.  However,  in  the  coherent  case  for  certain  angular  separations,  the  above 
inequalities  can  reverse  their  signs  and  in  that  case  the  above  threshold  expression  is 
meaningless.  This  is  because  the  relative  magnitudes  of  eigenvalues  can  be  different 
depending  upon  the  actual  angular  separation  and  consequently  the  additional  term 
B  (cj)  in  the  general  bias  expression  in  (51)  can  change  its  sign  so  as  to  reverse  the  ine¬ 
qualities  in  (8)  -  (10).  Since  B(u)  is  independent  of  N  and  a2,  from  (56.a),  for  large 
N ,  a  3  —  Q  ( ujm )  >  0  and  in  that  case  a  solution  always  exists. 

III.  Conclusions 

A  three-source  scene,  under  uncorrelated  and  coherent  situations,  is  separately 
analyzed  to  evaluate  the  sensitivity  and  robustness  of  MUSIC-type  high  resolution 
estimators  in  resolving  closely  spaced  sources.  When  the  null  spectrum  estimator  is 
used  to  locate  the  directions-of-arrival  of  incoming  signals,  for  a  fixed  number  of  sam¬ 
ples,  a  threshold  in  terms  of  SNR  exists  below  which  the  nulls  corresponding  to  the 
true  arrival  angles  are  no  longer  separately  identifiable.  Clearly,  in  a  two-source  scene 
the  corresponding  sources  are  separately  identifiable  if  the  bias  at  their  middle  angle 
is  larger  than  the  maximum  of  that  at  either  of  the  two  arrival  angles.  This  definition 
is  extended  to  three-source  scenes  here,  by  considering  pairs  of  signals  at  a  time  and 
choosing  the  maximum  of  the  respective  threshold  values  to  be  the  desired  SNR 
required  to  resolve  the  three  source  situation.  This  is  made  possible  by  first  obtaining 
parametric  expressions  for  the  eigenparameters  in  three  source  scenarios  and  using 
them  in  the  appropriate  bias  expressions.  The  parametric  expressions  for  eigenparam¬ 
eters  in  a  three-uncorrelated  and  coherent  situations  (the  two  outermost  signals  are 
perfectly  coherent  and  the  center  signal  is  uncorrelated  with  the  coherent  group) 
derived  here  is  believed  to  be  new. 

The  threshold  SNR  in  three  source  scenes  so  obtained  are  compared  with 
corresponding  results  in  two  source  situations.  These  results  indicate  that  in  the 
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uncorrelated  case,  the  presence  of  a  third  source  increases  the  SNR  required  to 
resolve  the  original  two  sources  approximately  by  a  factor  of  two.  It  may  be  remarked 
that  the  smaller  increment  in  SNR  in  the  ’symmetric’  coherent  case  can  be  attributed 
to  the  larger  angular  separation  (4u>d  instead  oi2ud)  between  the  coherent  signals  in 
the  three  source  case.  In  turn,  these  results  also  speak  for  the  sensitivity  of  MUSIC- 
type  eigenstructure  based  algorithms  when  an  extra  source  is  introduced  into  the 
scene. 

In  this  paper  we  have  only  analyzed  the  symmetric  coherent  case.  The  other  pos- 
sibility,  where  an  uncorrelated  signal  is  introduced  into  a  two  coherent  signal  scene 
such  that  the  angular  separation  between  the  uncorrelated  and  any  one  of  the 
coherent  signal  is  the  same  as  that  between  the  coherent  signals  deserves  further 
study.  This  ’more  natural’  case  seems  to  be  quite  involved  because  of  the  absence  of 
any  symmetry  that  is  present  in  the  coherent  case  anaylzed  in  this  paper. 
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Table  1 


Resolution  threshold  and  probability  of  resolution  vs.  angular  separation  for 
equipowered  sources  in  an  uncorrelated  scene,  (number  of  sensors  =  7,  number  of 
snapshots  =  100,  number  of  simulations  =  100).  Probability  of  Resolution  A  Total 
number  of  successes  in  100  simulations/100. 


j  Angles  of  arrival 

Angular 

Three  source 

Two  source 

Separation 

scence 

scence* 

Qi 

02 

e3 

2cod 

SNR(dB) 

Prob 

SNR(dB) 

Prob. 

23 

0.26 

9 

0.28 

34° 

o 

O 

45.33° 

0.1979 

24 

0.37 

10 

0.36 

25 

0.43 

11 

0.47 

26 

0.51 

12 

0.60 

14 

0.27 

3 

0.31 

60° 

66° 

71.73° 

022930 

15 

0.33 

4 

0.36 

16 

0.43 

5 

0.48 

17 

0.53 

6 

0.61 

11 

0.23 

2 

0.27 

127° 

135° 

144.33° 

0.3308 

12 

0.30 

3 

0.41 

13 

0.41 

4 

0.53 

14 

0.50 

*  In  the  case  of  two  source  scene,  the  first  two  angles  of  arrival  are  used  in  actual 
simulation. 


-28- 


0.15  025  035 

ANGULAR  SEPARATION 

Fig.  1.  Uncorrelated  source  scene:  Resolution  threshold  vs  angular  separation  for 
three  equipowered  signals  as  well  as  two  equipowered  signals.  A  seven  element  array 
is  used  to  receive  the  signals  in  both  cases.  Each  simulation  point  represents  one  hun¬ 
dred  trials  with  probability  of  success  for  resolution  ranging  from  0.33  to  0.5. 
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ANGULAR  SEPARATION  ( =  %T  omega  sub  d  %) 

Fig.  2.  Uncorrelated  source  scene:  Resolution  threshold  vs  angular  separation  for 
three  equipowered  signals  as  well  as  two  equipowered  signals.  A  ten  element  array  is 
used  to  receive  the  signals  in  both  cases.  Each  simulation  point  represents  one  hun¬ 
dred  trials  with  probability  of  success  for  resolution  ranging  from  0.33  to  0.5. 
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Table  2 


Resolution  threshold  and  probability  of  resolution  vs.  angular  separation  for 
equipowered  sources  in  a  symmetric  coherent  scene,  (number  of  sensors  =15, 
number  of  snapshots  =  100,  number  of  simulations  =  100).  Probability  of 
Resolution  A  Total  number  of  successes  in  100  simulations/100. 


Angles  of  arrival 

Angular 

Three  source 

Two  source 

separation 

scence 

scence* 

0i 

e2 

e3 

203^ 

SNR(dB) 

Prob. 

SNR(dB) 

Prob. 

28 

0.31 

23 

017 

155° 

158° 

161.45° 

0.0656 

29 

0.34 

24 

0.35 

30 

0.41 

25 

0.67 

31 

0.54 

19 

0.27 

12 

0.16 

110° 

112° 

114.03° 

01024 

20 

0.37 

13 

0.36 

2L 

0.48 

14 

0.59 

22 

0.64 

12 

0.25 

4 

0.23 

55° 

58° 

60.90° 

01372 

13 

0.37 

5 

0.38 

14 

0.47 

6 

0.63 

15 

0.61 

*  In  the  case  of  two  source  scene,  the  first  two  angles  of  arrival  are  used  in  actual 
simulation. 
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ANGULAR  SEPARATION  ( =  %T  omega  sub  d%) 

Fig.  3.  Coherent  source  scene:  Resolution  threshold  vs  angular  separation  for  three 
equipowered  symmetric  (two  coherent  and  one  uncorrelated)  signals  as  well  as  two 
equipowered  signals.  A  fifteen  element  array  is  used  to  receive  the  signals  in  both 
cases.  Each  simulation  point  represents  one  hundred  trials  with  probability  of  success 
for  resolution  ranging  from  0.33  to  0.5. 
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Appendix  A 


Consider  an  uncorrelated-source  scene  where  at  most  two  signals  may  be  com¬ 
pletely  coherent.  Thus  either  all  signals  are  uncorrelated  with  each  other  or  two  are 
perfectly  coherent,  while  the  rest  is  uncorrelated  with  each  other  and  with  the 
coherent  group.  In  both  situations  we  will  show  that  in  the  case  of  a  uniform  array,  if 
0V  02,  •  *  * ,  PM  denotes  one  set  of  eigenvectors  of  the  smoothed  covariance  matrix  R, 
then  their  inverted  and  complex  conjugated  counterparts  7.  =  J 0*,  i  =  1,  2,  •  •  • ,  M 
also  form  a  new  set  of  eigenvectors  of  R.  Here  by  definition 


Proof 

Let 


O  1 

1 

1  o 

. 


The  smoothed  covariance  matrix  has  the  form 

R  =  ARu  Af  +  ol 

where 


(A.1) 

(1.2) 


(A.3) 


1 

1  •  • 

1 

*1 

*2  •• 

*12 

s2  • 

*1 

ri 

0-y  •  ' 

tm 

6k 

(A.4) 


with  Sk  =  e  1  TW\  k  =  1,  2,  •  •  • ,  K  and  Ru  represents  the  smoothed  source  covari¬ 
ance  matrix.  Clearly,  in  the  uncorrelated  case 

Ru  =  diag  [P v P 2,  —,P£)  (A-5) 

and  in  the  latter  case  consisting  of  two  coherent  and  uncorrelated  signals  -  say  -  ( the 
first  and  the  second) 

Pl  pt  0  0  0  ' 

/>;  P2  0  0  0 

0  0  P3  0  0  (A.6) 

0  0  0-0 

0  0  0  0  ?^ 


where 

<2'‘,  •,^1] 

and  hence  (A.8)  simplifies  into 


(A.7) 

(A.8) 

(A.9) 
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GAGf  =  A($A/-1Rj$'(Af-1))At  +  a\ 


=  AR„  Af  +  a\  =  R 


(A.  10) 


since  1Rj$-(M_1)  =  R 


=  Rm  for  either  forms  of  Ru  given  by  (A.5)  or  (A.6). 


Interestingly,  if  some  of  the  entries  in  A  are  distinct,  say  the  first  K ,  then  since 


R  =  GAG*  =  BAB1 


with  V  =  B*G,  we  have 


VA  =  AV 


A ivU  =Xivij- 


Since  V  is  unitary  this  simplifies  into 


(A- 11) 


%  ~  P,  e  ‘  »  /  *  1,  2,  • 


(A- 12) 


7*T+2»  ***  “  [&K  +  1’  &K  +2> 


where  Vx  is  an  (Af  -K)  x(M  -K)  unitary  matrix. 


(A.  13) 
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Appendix  B 

In  this  Appendix,  we  evaluate  the  bias  expression  for  symmetric  coherent  signal 
scene  described  in  HC.  For  K  =  3,  the  first  term  of  the  general  bias  expression  in 
(48)  can  be  written  as 


3  f  3  lkldi 

«i  -  E  E  “ — — 
{;!  *  -A‘> 


(|#a(w)|2-  |#a(u,)|2) 


M  Lkkii 

+  S  — : 
‘■4  -\Y 


( I  A*a(“)  I  2  -  l#a(w)|2) 


(B.l) 


Then,  by  using  (A  12),  for  i ,  k  <  3  in  (50)  we  have  =  fiikk  and  for  /  <  3  and  k  >  4 

*  A  +  A,Rt  A)  -  2  A  • 

With  these  identities,  the  first  term  in  (B.l)  can  be  shown  to  be  zero  and  the  second 
term  can  be  expressed  as 


E 


11  2(A.  -<r7 


((A/  -3)  |#a(w)|  2-(2(w>).  (B.2) 


Now,  consider  the  second  term  of  (48)  by  partitioning  it  as 


3  f  3  3 


%**(")-  E  E  E 


,-i  A-1/-1  (A.  -A.) (A. 

kfii  tfi  ‘  k  1 


»*(«)£*  ^/a(w) 


3  Af  iM« 

+  E  E  — r 

k~i,m 4  (\  -\)(\ 


aW*  #a(u/)  +  S  E 


*.4 /-i  (A.  -A*)  (A,. 


•f(4  #/a(cj) 
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(B.3) 


M  M  ^ktii  +  -  -  + 

+  E  E  - aV)4W 

*;4'-4  -\)(\  -\) 


By  applying  the  orthogonality  property  between  Pk,  4<k  <M  and  a^- ),  1  < i  <  3  in 
(50)  and  with  the  help  of  (All)  -  (A  12),  it  is  easy  to  show  that  the  last  three  terms 
of  (B.3)  are  zeros.  Also,  from  the  structural  properties  of  and  Rb ,  we  obtain 

R*  =  +  /  2 P  M  sin2(Af  -  l)ud  ( uxuj  -  u3u| ) .  (B.4) 


Using  this  together  with  (50),  we  can  evaluate  the  first  term  of  (B.3).  For  /,  k,  l  <3 
(i,k,  l  are  all  distinct), 


4fM,  =  ffl  Rf  0  i>y  ’0,  +  0l  R‘  0,  0  R"  0,  -  20l  R7  0. 1 \ I  Rf  0, 


=  (  0l  R6  0t  0j  ( “1U3  “  )  0i  ]  -201^00*^  0t  .  (B.5) 

-  j  fa 

Here  we  have  made  use  of  the  fact  that  \  =  0.  e  and 

0,  =  2 A.  S„  -  0}r!  0  =  -  fit/  0,  for  i  fl  . 

Making  use  of  the  explicit  forms  of  0i ,  i  =  1,  2,  3  in  (45)  and  (47),  it  can  be  shown  that 

AW  -  u3uM-  -  0  »  **1.2,3  (B.6) 

and 

ty  % = ty  0 3=o.  IB.?) 

With  (B.6),  finally  (B.5)  reduces  to 

f m  =  -j0y0,0’R/0,.  (b.s) 

This,  together  with  (B.7)  and  (B.S),  simplifies  (B.3)  into 

333  rw,  .  . t 

£  £  — — : - >W,W 

(\  -\> 

kfil 
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(B.10) 


Re(j5t,V  fa \u)03fila(u))  . 


Using  (45)  and  (47),  since 

03V  02  =  j  2P M sm2(M -lyfa^l'^ 
and 

0*  R;  £i  =  -j  2PM  sin  2 (M  - 1)  jfl£  ux  u  J  pl , 
(B.IO)  can  be  expressed  as 

4P2M2sin22(M  -l)ud 


B(u)  = 


—  |  p\  ux  |  2  Re  ( a\u)fala(u) )  .  (B.  1 1) 


(A2-A1)(A2-A3) 

Finally,  with  (B.2)  and  (B.ll)  we  have  the  desired  bias  expression  as 


fj(u)  = 


Vx  +  V2 


1 


N  +  °V> 


3 

=  E 


\  2 
A-  a 


«•-!  2N  (A.  -  <x2)2 


((M  -3)  |#a<uO|2- £(*)}  +  ~  +  °(p-) 


In  a  similar  manner,  the  variance  expression  in  (49)  can  be  partitioned  as 


(B.12) 


Var{Q(u))  -~Y,  2 
i-1 1 


r  3  3 


3  M 


E  E  7,;/w  +  E  E  Qijkl 

k  =  \  /  » 1  &  *  1  /  =*4 


M3  M  M 

+  E  E  Qijkl  +  E  E  Qijkl 


Xc  -4  / -1 

ifi 


k  ”4  /  ■  4 


+  °(^ 


(B.13) 


where 


Re(  $kiji  at(w)3y  ^/a (w)  +  f^/f  aV)  A  ^a(u))af(u)4  /?/  a(w)  ] 

%ki  ~  _ 

(A,  -  A,) (A.  -  A,) 

From  the  symmetry  of  the  subscripts,  the  first  term  in  (B.13)  can  be  shown  to  be  zero. 
Moreover,  since  (Rb  )0k  =  0  for  i  <  3  and  k  >  4,  using  the  definition  of  f ijkl  in 
(50)  the  second  and  third  terms  in  the  above  expression  also  turns  out  to  be  zero. 

Now,  we  evaluate  the  fourth  term  in  (B.13).  For  <3  and  k ,  /  <  4,  we  have 

4*W  “  ' h  +  ^R‘  A)  -  2\  <r\  S,J  (B.14) 

and 

4  V  -  2\  %  0*'  ^  %  +  4*  K4  1,)  -  2va  A,.  „2  ,  (B.15) 

where  we  have  made  use  of  the  identities  \  =  J/?*  =  for  i  <  3  that  follow  from 
(45)  and  (47)  and 


-#tf(orR*)(S  Vffy-Jv,,  for  f>4. 

/>“  4 

Here  vw  is  an  element  of  Vr  With  (B.14)  and  (B.15)  in  (B.13),  the  variance  can  be 
simplified  as  [13] 

A  2  3  3  M  M  1 

Var(Q{u))  =  E  2  E  E  + 

/  “  1  /  ■  1  ^»4  /  »  4 


20(w)  3  V  -t  ,2  1 

-  — —  E  - 77  I  #«(■»)  I  +  O(-y) .  (B.16) 

^  ‘  -1  (A-  -  a2)2  N 

Notice  that  this  expression  has  structually  the  same  form  as  that  for  the  three 
uncorrelated-source  scene  in  (7). 
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Appendix  B 

A  New  Spectrum  Extension  Method  That 
Maximizes  the  Two-Step  Prediction  Error  - 
Generalization  of  Maximum  Entropy  Method* 


Abstract 


Given  (n+1)  consecutive  autocorrelations  of  a  stationary  discrete-time 
stochastic  process,  one  interesting  question  is  how  to  extend  this  finite  sequence  so 
that  the  power  spectral  density  associated  with  the  resulting  infinite  sequence  of 
correlations  is  nonnegative  everywhere .  It  is  well  known  that  when  the 
Hermitian  Toeplitz  matrix  generated  from  the  given  correlations  is  positive- 
definite,  the  problem  has  an  infinite  number  of  solutions  [1]  and  the  particular 
solution  that  maximizes  entropy  results  in  a  stable  all-pole  model  of  order  n. 
Since  maximization  of  entropy  is  equivalent  to  maximization  of  the  minimum 
mean  square  error  associated  with  one-step  predictors  [2],  in  this  paper  the 
problem  of  obtaining  admissible  extensions  that  maximize  the  minimum  mean 
square  error  associated  with  k-step  (k  <  n)  predictors,  that  are  compatible  with  the 
given  correlations,  is  studied.  It  is  shown  here  that  the  resulting  spectrum 
corresponds  to  that  of  a  stable  ARMA  (n,  k-1)  process.  The  details  of  this  true 
generalization  of  the  maximum  entropy  extension  are  worked  out  here  for  a  two- 
step  predictor  along  with  several  other  interesting  conclusions. 


I.  Introduction 

An  interesting  problem  in  the  study  of  autocorrelation  forms  and  their 
associated  power  spectral  densities  is  that  of  estimating  the  spectrum  from  a 
finite  extent  of  its  correlation  function.  Known  as  the  trigonometric  moment 
problem  in  the  discrete  case,  it  has  been  the  subject  of  extensive  study  for  a  long 
period  [1  -  6].  In  view  of  the  considerable  mathematical  interest  as  well  as  the 
practical  significance  of  the  moment  problem  in  interpolation  theory,  system 
identification,  power  gain  approximation  theory  and  spectrum  estimation,  it  is 
appropriate  to  review  the  problem  briefly  in  the  present  context.  Towards  this,  let 

*  This  work  was  supported  by  the  Office  of  Naval  Research  under  contract  N-00014-89-J-1512. 
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x(nT)  represent  a  discrete-time,  zero  mean,  wide  sense  stationary  stochastic 
process  with  autocorrelation  function, 

ifc  =  E [x(nT)  x*((n  +  k)  T)]  =  r^,  k  =  0,l,2,-..  -  .  (1) 

As  is  well  known,  the  power  spectral  density  S(0)  of  this  stationary  process  is 
given  by  the  discrete-time  Fourier  transform  of  its  autocorrelation  sequence  [1], 
i.e, 

m  =  £  rk  eJ“ .  (2) 

k»-« 


Moreover,  S(0)  £  0  and 


^  fs(0)  e-^d©,  Ikl  >0  . 


Clearly  for  processes  with  finite  power,  we  have 


(3) 


1_ 

2k 


S(0)d0=  r0  =  E[(x(t)|2; 


<°°, 


(4) 


i.e,  S(0)  is  integrable  (belongs  to  over  -k  £  0  £  k).  The  non-negativity  property  of 
the  power  spectral  density  can  be  characterized  in  terms  of  certain  Toeplitz 
matrices  generated  from  its  correlations  [2].  Let  Tn  denote  the  Hermitian  Toeplitz 
matrix  generated  from  r0,  rj,  r2, .  . .  rn,  and  An  its  determinant.  Thus, 


* 


♦ 


Ln-1 


♦ 


(5) 
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and(1) 


Then  [2,  6] 


An  =  det  Tn  # 


S(6)£  0  A„  £  0  ,  n  =  0,l,2,-~  , 


(6) 


i.e.,  the  nonnegativity  property  of  the  power  spectral  density  function  is  equivalent 
to  the  nonnegative  definiteness  of  all  Hermitian  Toeplitz  matrices  generated  from 
its  correlations. 

Further,  assume  that  the  process  also  satisfies  the  Paley- Wiener  criterion 
(causality  criterion)^) , 


H  = 


2ft 


£ 


In  S(0)d0>-< 


(7) 


i.e. ,  the  entropy  H  of  this  process  is  finite.  With  this  additional  constraint,  it  can  be 
shown  that  [7  -  8]  the  nonnegative  property  of  the  power  spectral  density  implies 
positive-definiteness  for  all  Tn  in  (5),  i.e. ,  An  >  0,  n  =  0, 1,  2, . . . , 00 . 

The  integrability  condition  in  (4)  together  with  the  Paley- Wiener  criterion 
permits  the  factorization  of  the  power  spectral  density  in  terms  of  a  unique 
function  with  certain  interesting  properties.  In  fact,  in  that  case,  there  exists  a 
unique  function  [7  -  9] 

00 

B(z) »  £  bkZk,  bo>0  (8) 

k-0 

that  is  analytic  together  with  its  inverse  in  |  z  |  <  1 ,  such  that 


(!)  Here  onwards  a  (or  A),  a  and  A  stand  for  scalar,  vector  and  matrix  in  that  order. 
Similarly,  A*,  and  A^  represent  the  complex  conjugate,  transpose  and  complex  conjugate 
transpose  of  A,  respectively.  The  symbol  det  A  ■  I A I  is  used  to  denote  the  determinant  of  the 
matrix  A. 

(2)  From  (4),  the  inequality  H  <  +  «  is  automatically  satisfied  [8]. 
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and<3) 


X  l^kl  <°° 

k-0 

S(0)=lB(eje)l2>  a.e. 


(9) 


(10) 


This  minimum  phase  factor  B(z)  (free  of  zeros  and  poles  in  |z|  <  1)  is  known  as 
the  Wiener  factor  of  the  given  process  and  represents  a  causal,  digital  filter  with 
square  summable  impulse  response.  Its  physical  meaning  is  not  difficult  to 
grasp.  When  driven  by  an  appropriate  stationary  white  noise  source  of  unit 
spectral  density,  this  filter  regenerates  the  given  stochastic  process  x(nT),  entirely 
from  the  past  samples  of  the  input  white  noise  process  (see  Fig.  1). 


w(nT) 


x(nT)  =  £  bfc  w(  (n  -  k)  T  ) 
k  =  0 


Fig.  1.  Wiener  Filter  of  a  Stationary  Process 


We  are  in  a  position  to  state  the  spectrum  extension  problem:  Given  (n  +  1) 
autocorrelations  r0,  rj.,  t2, . . .  rn,  from  a  stationary  stochastic  process  that  satisfies 
(4)  and  (7),  determine  all  solutions  for  the  power  spectral  density  S(0)  that  are 
compatible  with  the  given  data;  i.e.,  such  a  solution  S(0)  should  satisfy 

S(0)  >  0  (11) 

and 


2k 


S(0)  e-*6  d0  =  iv, 


k  =  0, 1,  2,---  n  , 


(12) 


B(e*9)  *  lim  B(re*9) .  In  addition  to  aesthetic  reasons,  the  use  of  the  variable  z  for  the  delay 

r  -» 1  -0 

operator  (in  contrast  to  the  usual  z~l),  translates  all  stability  arguments  to  be  carried  out  in  the 
compact  region  |  z  |  S  1.  Notice  that  r±,  k  =  0  ->  «>  real  guarantees  b*,  k  =  0  «» to  be  real;  i.e.,  the 

Wiener  factor  for  real  processes  are  real. 
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in  addition  to  satisfying  (4)  and  (7). 

It  is  well  known  that  a  necessary  and  sufficient  condition  for  the  existence 
of  such  a  solution  is  the  nonnegativity  property  of  the  Toeplitz  matrix  T„ 
generated  by  the  given  (n  +  1)  correlations  [2,  3].  Moreover,  when  Tn  is  positive 
definite  (such  is  the  case  here),  the  above  spectrum  extension  problem  has  an 
infinite  number  of  solutions(4)  [3,  4].  Youla  has  parametrized  this  entire  class  of 
solutions  for  the  spectral  extension  problem  in  a  network  theoretic  setting  in 
terms  of  bounded  real(5)  (b.r)  functions  [10].  These  solutions  also  follow  from 
Schur's  theory  on  b.  r.  functions  [3].  Before  examining  them,  to  make  further 
progress  in  the  present  context,  it  is  necessary  to  gain  some  additional  insight  into 
the  concepts  such  as  maximization  of  entropy,  prediction  errors,  and  the  Wiener 
factor. 

Section  II  examines  the  maximum  entropy  extension  and  discusses  its 
significance  in  terms  of  possessing  a  Wiener  factor  that  maximizes  the  minimum 
mean  square  error  associated  with  one-step  predictors.  Further,  Youla’s 
parametrization  of  the  class  of  all  extensions  that  are  compatible  with  the  given 
data,  and  their  dependence  on  the  maximum  entropy  solution  [10]  is  briefly 
reviewed. 

Given  r0,  rlf . .  .  rn,  the  Wiener  factor  that  maximizes  the  minimum  mean 
square  error  associated  with  a  k-step  predictor  is  shown  to  have  an  ARMA  (n,  k- 
1)  structure.  The  special  case  of  the  two-step  predictor  is  treated  in  Section  III 
and  using  an  induction  argument  the  general  case  is  proved  in  Appendix  A. 


When  ^  =  0,  the  above  extension  problem  has  a  unique  solution. 

(5)  p(z)  is  said  to  be  bounded  real  (b.r)  if  I  p(z)  I  SI  in  I  z  i  <1  and  is  real  for  real  z.  By  Cauchy's 

am 

inequalities  if  p(z)  =  £  pkzk  ,  then  I  pk  I  <  1  for  all  k.  Further,  if  I  p0  I  =  1  ,  then  p(z)  *  1  . 
k»0 
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Using  Youla’s  parametric  formulation,  sin  explicit  solution  for  the  two-step 
predictor  is  worked  out  and  presented  in  Section  III. 


II.  Maximum  Entropy  Solution  and  the  Class  of  All  Extensions 

A  geometric  interpretation  of  the  class  of  admissible  solutions  is  well 
known  [1,  4,  5].  Clearly,  under  the  hypothesis 


To 

rl 

*2 

"•  rn 

^  = 

* 

rl 

r0 

rl 

rn-l 

* 

rn 

Cl 

•  *  ‘ 

ro 

any  admissible  extension  r^,  k  »  n  +  1  -» <*>,  should  satisfy 


Ak  >  0,  k  =  n  +  l-»' 


Consequently,  at  the  first  step  rn+1  =  x  must  be  chosen  so  that 


An+i(x)  = 


fo  rl 


ri  r0  *i 


*  *  * 
x  r„  r  , 
n  n-1 


rn  x 


>  0 


(13) 


(14) 


(15) 


(16) 


Using  well  known  matrix  identities*6*  7> ,  An+i*x)  can  be  expanded  as 

4,-1  l*-5nl2. 

where 

^“YfT^Yb  (17) 

Yf  =  [rl>  r2 . rn]T  . 

and 

Yb^Wn-i . rJT. 

In  turn,  (16)  implies  that,  the  set  of  values  x  that  satisfy  (15)  is  the  interior  of  a 
circle  with  radius 

^  =  -^->0  (18) 

Vi 

and  center  given  by  (17). 

The  general  rule  is  now  conceptually  clear:  Having  selected  rn+i  =  x 
consistent  with  (15),  rn+2  can  be  chosen  from  the  interior  of  a  circle  with  radius 
Xn+i  =  An+]/An  and  center  2;n+1  and  so  on.  Notice  that  except  for  Xn  and  which 
are  uniquely  determined  by  the  given  data,  all  successive  and  ^  depend  on  the 
particular  rule  used  to  pick  r^  from  the  interior  of  the  respective  circles  k  =  n+2  -> 


(6)  Let  A  be  an  n  x  n  matrix  and  Anw,  A^,  AgW,  A^  denote  the  (n  -  1)  x  (n  -  1)  minors  formed 

from  consecutive  rows  and  consecutive  columns  in  the  northwest,  northeast,  southwest  and 
southeast  comers.  Further  let  Ac  denote  the  central  (n  -  2)  x  (n  -  2)  minor  of  A.  Then  from  a 

special  case  of  an  identity  due  to  Jacobi  [11], 


Ag  I A I  —  Anw  Agg  Ane  A, 


sw 


(7) 


A  B 
C  D 


|  A  |  |  D  -  CA-lD  | 
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oo.  One  particular  choice,  however,  has  several  interesting  properties.  To  observe 


this,  first  notice  that^  for  any  k, 


a-^(n2-|4<1>|2 


*k+l 


) 


(19) 


which  gives  the  useful  identity 


^k+1  = 


Ak+1  Ak 


Ak  Ak-1 


1  - 


Al) 

^k+1 


Ak 


2\ 


^  Xt  . 


(20) 


Here,  by  definition  A^l  denotes  the  minor  of  A^+1  obtained  by  striking  out  its  first 


column  and  last  row.  From  (20),  the  sequence  X^  of  positive  numbers  for  k  =  n 
oo  is  always  monotone  nonincreasing  and  has  a  limit  given  by  [7] 

Ak 


lim  Xu  = 

Atj  0 


=  K  >0 


(21) 


with  b0  as  in  (8)  representing  the  constant  term  in  the  Wiener  factor.  Since  X^  <  X„ 
for  k  =  n  ->  oo,  this  further  implies  that 


with  equality  if  and  only  if 


Or,  equivalently,  in  that  case 


x  =  0  ,  k  =  n->oo. 


.  an  , 

^  =  — ,k  =  n 


(22) 


(23) 


From  (22),  (23)  and  (16),  b0  assumes  its  maximum  possible  value 

(bo)r 

if  and  only  if 


-'o-'max 


An_ 

^-1 


(24) 
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(25) 


x  =  rk+1  =4.  k  =  n -> °<= , 

i.e.,  only  when  all  extensions  are  at  the  centers  of  the  respective  admissible  circles. 
2 

Since  b0  represents  the  minimum  mean  square  error  associated  with  a  one- 
step  predictor  that  makes  use  of  the  entire  infinite  past(8)  [2],  the  above  completion 
rule  suggests  that  the  unique  extension  so  obtained  by  identifying  each  r^,  k  =  n  + 
1  -»  oo  with  the  centers  of  the  admissible  circles  also  maximizes  the  minimum 
mean  square  error  of  a  one-step  predictor  of  the  given  process,  i.e.,  of  all  the  one- 
step  predictors  that  can  be  generated  from  each  admissible  completion,  the  one 
above  has  the  maximum  value  for  the  minimum  mean  square  prediction  error 
and  in  that  sense  it  is  maximally  robust.  Interestingly,  this  is  also  the  maximum 
entropy  extension  originally  formulated  by  Burg  [12,  13].  To  see  that  the  above 
extension  also  possesses  maximum  entropy  among  all  possible  extensions,  it  is 
enough  to  relate  the  entropy  H  of  the  process  and  the  prediction  error  b0. 

Towards  this  purpose  let  d^,  ikl  s  0,  1,  2,  .  .  .  represent  the  Fourier 
coefficients  of  In  S(0).  Thus,(9) 

dk  =  i_  £  ^  S (0)  e_jk9  d0  (26) 

and  hence  [4] 

oo 

In  S(0)  =  X  4^®  (27) 

k  =  -«> 
or 

^  Using  least  mean  square  theory  it  follows  that  the  best  one  step  linear  predictor 

n-1 

A 

x(T) »  £  wk  x{-kT)  that  makes  use  of  n  past  samples  results  in  a  mean  square  error 

k-0 

1 4„_!  Clearly  from  (21),  lim  8„  =  bo  . 

(9)  The  Paley-Wiener  criterion  (7)  guarantees  that  lnS(Q)  is  well  defined  almost  everywhere. 
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S(6)  =  exp( d0)  exp  X 

k«l 


where 


exp  I  dke*e  =L(eie)L*  (e^9)  (28) 

k-1 


L(z)  a  exp  (do/2)  exp  X  dkzK  A  exp  (do/2)  X  ak  zK . 
~  \  k  =>  1  —  k  =  0 


Clearly  a0  =  1.  On  comparing  (10)  and  (28),  we  can  associate(10)  L(z)  with  the 
Wiener  factor  B(z),  and  in  that  case  by  equating  (8)  and  (29),  we  obtain 

bk  =  ctfe  exp  (do/2) ,  I k|  £  0  (30) 


In  particular,  from  (26)  [2,  5,  6] 


bf  =  exp(do)  =  exp 


if 

2nJ„ 


In  S(0)d0  a*  exp(H) 


This  one-to-one  relationship  between  entropy  H  and  the  one-step  prediction 

2 

error  allows  one  to  conclude  that  the  above  extension  method  that  maximizes  b0  is 
also  the  maximum  entropy  extension.  The  maximum  entropy  solution  plays  a 
basic  role  in  parametrizing  all  other  admissible  extensions.  In  a  network 
theoretic  setting  involving  positive-real  (p.r.)  functions,  Youla  has  shown  that,  in 
the  case  of  real  correlations,  every  admissible  solution  S(0)  can  be  represented  as 


S(e)=|Bp(ej0)f 


where,  the  Wiener  factor  B0(z)  is  given  by 


«  M  Hz) 

b'(z,=d M 


Here,  Hz)  is  a  b.r.  function,  analytic  together  with  its  inverse  in  |z|  <1  (minimum 


(10)  Notice  that  L(z)  is  also  analytic  together  with  its  inverse  in  Izl  Si. 
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phase  factor),  that  is  associated  with  the  spectral  factorization 

l-|p(eje)|2*|r(eje)|2 

or,  more  compactly,  with  the  factorization 


where,  by  definition 


1  -  p(z)  p*(z)  *  Hz)r.(z) 


p*(z)  =  p(l/z) 


and  p(z)  is  an  arbitrary  bounded  real  function^.  Further, 

Dn(z)Apn(z)-zp(z)Pn(z)  (35) 

and  Pn(z)  is  the  unique  degree  n  polynomial  generated  by  the  Levinsion  algorithm 
from  the  given  correlations  r0,  rj,  . . .  rn,  through  the  recursion  [14]: 

VT-s£  Pk(z)  =  Pk4(z)  -  skzk Pk4(l/z) ,  k  =  1 , 2,  •  •  •  (36) 

Here  sk,  k  =  1 ,  2, . . .  are  the  reflection  coefficients  given  by 


sk  =  (-1) 


k-l  k 


with  Ak  as  defined  in  (19M20).  This  recursion  is  carried  out  under  the 


initialization 


Po(z)  = 


Interestingly,  Pn(z)  can  be  compactly  written  as(11) 


(11)  The  Levinson  recursion  (36)-(37)  follows  from  (38)  through  an  easy  determinant 
expansion^. 
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ro  rl 


Pn(2)  = 


Vv^7 


*1  To 


rn-l 

„n  ,n-l 


rn-l 


ro  1*1 


A  ao  +  ajz  +  a2Z2 +  •••  anZn  (38) 


and  using  Rouche's  theorem  [15]  (or  otherwise)  on  (36)  and  (35)  repeatedly, 
through  an  induction  argument  it  is  easy  to  show  that  Pn(z)  and  hence  Dn(z)  are 
strict  Hurwitz  Polynomials.(12)  Moreover,  by  definition  Pn(z)  represents  the 
polynomial  reciprocal  to  Pn(z),  i.e. , 

Pn(z)  =  znPn(l/z). 


Thus,  p(z)  b.  r.  implies  that  Bp(z)  is  analytic  together  with  its  inverse  in 
|  z  |  <  1  with  square  summable  Fourier  coefficients.  Clearly,  p(z)  parametrizes  S(0) 
and  moreover  all  these  extensions  satisfy  the  correlation  matching  property 
involving  the  first  (n+1)  coefficients.  In  Youla's  representation  [10],  all  such 
spectral  extensions  can  be  realized  as  the  real  part  (on  the  unit  circle)  of  the  input 
impedence  of  a  cascade  of  (n+1)  lossless,  equi-delay,  transmission  lines  that  has 
been  terminated  upon  an  arbitrary  passive  load  W(z).  The  transmission  lines 
with  characteristic  impedences  R0,  Ri, .  . .  ,Rn,  are  generated  from  the  given  real 
correlations  such  that,  the  junction  (mismatch)  reflection  coefficients  S]c=(Rk~Rk- 
1)/(R]C+R]C.1) ,  k  =  1, 2,  •  •  •  n  are  given  by  (37)  with  R0=r0.  In  this  representation  p(z) 
represents  the  reflection  coefficient  of  the  load  W(z)  at  the  far  end  normalized  to 
the  characteristic  impedence  Rn  of  the  last  line,  i.e.. 


A  Hurwitz  Polynomial  has  all  its  zeros  in  Izl  >  1.  A  strict  Hurwitz  polynomial  has  ail  its 

zeros  in  Izl  >1.  Note  that  r*,  k  =  0  -*  n  real  implies  that  the  coefficients  a*,  k  =  0  -»  n  in  (38) 
are  also  real. 
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(39) 


p(z)  = 


W(z)-Rn 


W(z)  +  Rn 

To  obtain  the  maximum  entropy  solution  explicitly,  we  proceed  to  evaluate 
the  entropy  Hp  associated  with  the  general  solution(32)  and  maximize  it  with 
respect  to  the  free  parameter  p(z).  Since 


Hp 


=  -f 

2*J-„ 


In  S(0)  d0  =  In  bg  , 


from  (8)  and  (32) 


b*  =  Bp(0)  =  =  r2(0) 

J  P  D*(0)  P*(0)  Vi 


(40) 


and  using  (24),  (31) 


Hp  =  In  -p-  -  In  [l/  r2(0)]  =  Hme  -  In  [l/  r2(0)] 

-l 


p  Vi 


(41) 


Clearly  the  extension  that  maximizes  entropy  is  the  one  where  HO)  =  1 .  Since  Hz) 
is  also  b.r.,  notice  that  r2(0)  <  1  unless  Hz)  si,  and  consequently  HO)  =  1  implies^ 
p(z)  s  0.  In  Youla's  representation,  this  is  equivalent  to  terminating  the  last  line 
on  its  characteristic  impedence  Rn  (see  (39)).  Thus  from  (32)  -  (35)  the  maximum 
entropy  spectral  extension  has  the  form 

1 


S(0)  = 


P„(eje)l: 


(42) 


and  the  associated  Wiener  factor 

1 


Bme(z)  = 


Pn(z)  a0  +  axz  +  a2z2  + . . .  +  anzn 


(43) 


represents  a  stable  AutoRegressive  form  of  order  n  (i.e.,  AR(n)).  Alternatively,  van 
den  Bos  has  shown  that  [16]  the  standard  linear  prediction  method(8)  also  leads  to 
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the  same  set  of  Yule-Walker  equations  as  the  maximum  entropy  method  and 
hence  (43)  also  represents  the  best  linear  prediction  filter.  From  (32),  (33)  and  (35) 
the  polynomial  Pn(z)  that  characterizes  the  maximum  entropy  solution  plays  a 
key  role  in  all  other  extensions,  and  in  particular,  the  one  that  maximizes  the  two- 
step  prediction  error. 


m.  The  Wiener  Factor  that  Maximizes  the  Two-Step  Prediction  Error 

Given  the  correlations  r0,  r*,  r2, . .  .  rn,  of  all  the  admissible  completions 
given  by  (32)  and  (33),  the  problem  here  is  to  find  the  one  that  maximizes  the  k-step 
minimum  mean  square  prediction  error.  In  what  follows,  we  first  deal  with  the 
two-step  predictor  case.  It  is  shown  here  that  maximization  of  the  two-step 
prediction  error  results  in  an  ARMA  (n,  1)  process.  Through  a  constructive 
procedure,  the  existence  of  the  ARMA  (n,  1)  Wiener  factor  is  demonstrated  for  this 
case.  The  general  case  is  dealt  with  in  Appendix  A,  where  it  is  shown  that 
maximization  of  the  k-step  minimum  mean  square  error  results  in  an  ARMA  (n, 
k-1 )  process. 

Using  (29)  and  (30),  the  two-step  prediction  error  P2  is  given  by  [2], 


f71 

P2=  !  t>0l 2  +  1  bi  1 2  =  [1  +  |<Xil2]exp 

—  In  S(0)  d0 

2*  L 

Naturally,  maximization  of  P2  is  with  respect  to  the  unknown  autocorrelations 

1  f* 

rn+i>  rn  +2,  rn+3,  •  •  ■  and  using  the  relation  ax  =  dx  =  —  I  In  S(9)  e  J0  d6  ,  this  leads  to 

2*J-  it 


- 

dP2  _  1 

( 1  +  ai  |2)  +  ax  e->9  +  (Xi  e--*9 

“  2k 

4 

—  K 

S(0) 

ejk9d0 
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f * 


1_ 

2k 


1  +  dj  eJ0 


S(0) 


.  d0  =  0  ,  |  k  |  >  n . 


(45) 


-K 


Clearly,  (45)  implies  that  the  Fourier  series  expansion  for  the  real  periodic 
nonnegative  function  1 1  +  a*  ©i®  1 2/S(0)  truncates  after  the  nth  term  and  hence  it 
must  have  the  form 


1  +  a^9 


S(0) 


=  I  €**•- 

k*-n 


I  gkejke 

k  =  0 


(46) 


or 


m  =  ; 


2 

1  +  o^e^ 


=  |  B^ej0)  |2 


I  gke^ 

k»0 


(47) 


where 


Bg(z)  = 


A(z) 


G(z)  ’ 

A(z)  =  ( 1  +  z )  or 


'i-L.' 

* 

<*i 


(48) 

(49) 


and 


n 


G(z)  =  £  gkZk  . 
k*0 


(50) 


Thus,  at  least  formally  B2(z)  -  ARMA  (n,  1),  i.e.,  the  Wiener  factor  that 
maximizes  the  two  step  prediction  error,  if  it  exists,  is  of  the  type  ARMA  (n,  1).  To 
complete  the  argument,  we  must  demonstrate  the  existence  of  such  a  factor  that 
is  analytic  together  with  its  inverse  in  \z\  <1. 

Toward  this  purpose,  notice  that  in  the  case  of  real  correlations,  this 
specific  extension,  if  admissible,  should  follow  from  (32)  for  a  certain  choice  of  the 
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bounded-real  function  p(z)  and  in  that  case,  on  comparing  (48)  and  (33),  because  of 
degree  restrictions,  the  simplest  p(z)  must  have  the  form(13> 


p(z)  = 


a  +  bz 

For  (51)  to  be  bounded  real,  it  is  necessary  that  there  exist  no  poles  in  |  z|  <  1,  i.e., 


(51) 


>  1 , 


(52) 


and 


a  + 


beje 


£  1  <=*  (a±b)  2  1  , 


(53) 


(in  addition  to  a  and  b  being  real). 

Conversely,  when  (52)  and  (53)  is  true,  from  maximum  modulus  theorem 
[17],  analytidty  in  I zl  <1,  together  with  (53)  implies  boundedness  for  p(z)  in  Izl  S 
1.  In  the  present  situation,  the  existence  of  such  a  b.  r.  function  as  in  (51)  can  be 
verified  by  solving  for  a,  b  from  (34)  -  (35)  and  examining  whether  they  satisfy  the 
necessary  and  sufficient  conditions  (52)  and  (53).  In  that  case,  through  a  direct 
calculation,  the  degree  n  requirement  for  G(z)  yields, 


where  a,  P  satisfies 


r(z)  = 


a  +  Pz 
a  +  bz 


a2  +  p2  =  a2  +  b2-l 
aP  =  ab 


(54) 

(55) 


(56) 

(57) 


^■3)  It  is  shown  in  Appendix  B  that  p(z)  =  — 1Kf  ^  2  etc,  are  not  acceptable.  (Note  that  B0(z) 

£L  DZ  "f  CZ  r 

rational  implies  that  p(z)  must  be  rational.) 
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and 


Dn(z)  * 


I  8k  zk 

k»0 


a  +  bz 


(58) 


with 


gk  =  ak_!  b  +  ak  a  -  a,*_k+1 ,  k  =  0, 1, 2,  • .  •  n  . 


(59) 


Here  ak,  k  =  0, 1,  2,  .  .  .  n,  are  the  coefficients  of  Pn(z)  (see  (38)).  Moreover,  the 
Wiener  factor  B2(z)  yields  the  power  series  expansion  (in  I  zl  <  1), 


Ba(z)  = 


a  +  Pz 

2  Sk2l 
k=0 


a 

go 


A 

go  ’  So  / 


z  +  •  ■  •  4  b0  +  bjz  +  b2Z2  + 


(60) 


and  hence  from  (29)-(30) 


^  *>o  U  8o 


(61) 


On  the  other  hand,  from  (48)-(49)  and  (60) 


P  a 

*mS  or  J 


(62) 


(since  bfc ,  o^ ,  k  =  0  -* «>  are  real  in  the  case  of  real  correlations).  It  is  easy  to  show 
that  the  first  choice  ax  =  p/a  does  not  always  lead  to  a  bounded  real  solution  for 
p(z).  In  fact,  by  letting  ax  =  p/a  and  equating  this  to  (61),  we  obtain  gi  /  go  =  0, 
which  in  turn  implies  gi  =  0  since  go  =  aoa  is  a  finite  number,  i.e., 

g\  =  aob  +  axa  -  a,,  =  0 


or 


which  gives 


a  = 


,  2  2 
an  -  aob  -a0 


ai 


aia„ 
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a 

b 


2  2 

aoai 


From  (52),  for  p(z)  to  be  b.r.,  I  a  /  b  I  >  1 .  However  as  the  strict  Hurwitz  polynomial 
P2(z)  as  (z  +  2X2z  +  3)  =  2z2  +  7z  +  6  shows 

a  _  62  -  22  _  &  <r 
b  =  6x7  =  42  < 

and  hence,  from  (52),  as  (3/a  is  not  an  appropriate  choice.(14)  Turning  back  to  (62), 
this  leads  to  the  only  other  possibility 


a 


Equating  (61)  and  (63),  we  obtain 

p2-ct2  =  Si_ 
ap  ~  So' 


(63) 


(64) 


Using  (56),  (57)  and  (59)  and  after  some  algebra,  (64)  reduces  to  the  cubic  equation 
with  real  coefficients 


x3  +  px  +  q  =  0  , 


(65) 


where 


and 


p  as  -2 


1  + 


a 


2a 


2  \ 


<  0 


'o/ 


(66) 

(67) 


<14)  At  times,  the  above  choice  can  lead  to  admissible  solutions.  For  example,  in  the  case  of  the 

strict  Hurwitz  polynomial  P2(z)  =  5  -  2z  +  z2,  the  parameters  a,  b  so  obtained  generate  a  b.  r.  p(z) 
leading  to  an  admissible  Wiener  factor. 
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q  =  -^(l-i).  (68) 

Clearly,  the  Wiener  factor  B2(z)  in  (60)  that  maximizes  the  two-step  prediction 
error  represents  an  admissible  solution,  if  the  cubic  equation  (65)  has  at  least  one 
real  solution  with  magnitude  greater  than  unity  and  further  that  solution 
together  with  (54)  satisfies  (53).  To  examine  this,  notice  that  if  the  discriminant 


»-UMJ 


is  negative,  then  (65)  has  three  real  roots  and  if  D  >  0  it  has  one  real  root  and  two 
complex  roots  that  form  a  conjugate  pair  [18].  However,  as  shown  in  Appendix  B, 
the  above  discriminant  is  always  negative  and  the  corresponding  three  real  roots 
can  be  obtained  explicitly  by  making  use  of  Cardano's  formula.  In  that  case,  let 


R  =  sgn(q)l-£)1  =  -  sgn  (  ^)a/  f[i  +  !_  + A; 


Then,  the  three  roots  are  given  by  [18] 


xi  *  -  2  R  cos  ( <p  /  3  ) 

x2  =  -2Rcos((p/3  +  27i/3) 


X3  =  -2Rcos(<p/3  +  47c/3) 


where 


cos<p  =  ^ 


2  ,  1 
-  1  +  —  + 


.2  M  3/2 


It  is  easy  to  show  that  two  of  these  roots  always  have  magnitude  greater  than 
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unity.  This  follows  from  a  well-known  sufficiency  condition  due  to  Cohn(15), 
which  in  the  case  of  (65)  reduces  to  |p|>|q|  +  l,  for  two  of  its  roots  to  have 
magnitude  greater  than  unity.  To  verify  this  condition,  notice  that 


1=2 

a2  ) 

1+1  +  a‘. 

-2 

1 

•>2  2a„2  J 

ao 

b2 

-1 


>  0 


(73) 


and  hence  (65)  has  two  roots  with  magnitude  greater  than  unity  and  one  with 
magnitude  less  than  unity.  Moreover,  without  exception^16) ,  cos  cp  £  0  ,  which 
implies  |  <p  |  <  it  /  2  or 


cos 


>  cos 


(74) 


and 


1 

2 


<  -  cos 


(p  2  K  ) 

,3+tJ 


(75) 


Thus,  X2  and  X3  have  the  same  sign  that  is  opposite  to  the  sign  of  xj.  Clearly,  the 
sign  of  X2  and  X3  is  the  same  as  that  of  q.  The  structure  of  these  roots  is 
summarized  in  Fig.  2.  Further,  since  |  R|  >  V 2/3,  using  (74)  and  (75)  we  also  have 
2|R|>|x1|>V3|R|,  f3|Rj>|x2|>|Rj  and  |  x3 1  <|  R|  which  gives 


|  xi  |  >  |  x2 1  >  |  x3 1  •  (76) 

Thus,  using  the  Cohn  criterion,  X]^  and  x2  have  magnitudes  greater  than 
unity,  i.e .,  p(z)  obtained  using  any  one  of  these  roots  is  always  analytic  in  |  z|  £1  . 


(15)  por  a  polynomial  f(z)  =  ao  +  ajz  +  •  •  -  +  apZp  +  •  •  •  +  anzn  ,  ap  *  0  ,  if  |  ap  |  >  |  ao  j  + 1  ai  |  +  •  •  • 
+ 1  a^t  |  + 1  ap+1 1  + .  ■ .  + 1  a„  |,  then  f(z)  has  exactly  p  zeros  inside  the  unit  circle  [15]. 

(16)  |  b  |  =  |  ao  /  a„  j  >  1  ,  since  a0  /  a„  represents  the  product  of  the  n-roots  of  the  strict  Hurwitz 
polynomial  Pn(z). 
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For  bounded  reality,  it  remains  to  show  that  a,  b  so  obtained  from  these  roots  also 
satisfy 

(a  ±  b)2  >  1. 


In  what  follows,  we  will  first  demonstrate  this  condition  for  the  largest  rootf17>  Xj. 
In  that  case,  using  (66),  (70)  and  (74),  x*  gives  rise  to 


,<i> 


*i 


b|  £  V3|R||b|  =  a  2 


i+b2+-sl 

2aSj 


>  2 


Again,  without  loss  of  generality,  we  can  assume  a<:i)  and  b  are  of  the  same  sign  - 
say  positive  ( i.e. ,  x1  >  0  ,  b  >  1  )  -  and  hence, 

( a(1)  +  b  )2  >  9  >  1 . 

To  verify  the  second  part,  let 

\2 


(77) 


(  a(1)  -  b  )2  >  ^j~2 


1  +  b2  +  ~~z  -b 
2a;/ 


*  A 


(78) 


3A 


To  find  its  minimum,  notice  that  by  straightforward  algebra  —  =  0  yields 

ob 


u2  i 

bn  =  1  +  — . 


But,  the  second  derivative 


92A 

9b2 


n 


'  a-,2  ^ 


1  + 


l 


l+b2  +  -^ 


n* 


Since  |  x3  |  <  1,  there  can  atmost  be  two  solutions,  one  corresponding  to  and  the  other 
corresponding  to  *2-  The  solution  due  to  *2  is  discussed  later. 
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~2  ax 

is  always  positive  and  hence,  in  particular,  at  bo  =  1  +  — - ,  the  quantity  A  in  (78) 

24 


achieves  its  minimum  value  given  by 


or 


in]in=V2(2b„2)-b0  =  b0 

(a(1)-bo)2  2ii  =  bi>l. 


This  completes  the  proof.  The  uniqueness  for  this  choice  of  p(z)  is  proved  in 
Appendix  C. 

To  investigate  the  solution  due  to  X2,  notice  that  under  the  above 

/n\ 

assumptions  (  X!  >  0,  b  >  1  ),  X2  and  hence  &  =  x2  b  are  negative  and  as  a  result 

|  a(2)  -  b  |  >  1  is  automatically  satisfied.  For  bounded  reality  of  the  p(z)  so  obtained 
from  a(2)  and  b,  it  remains  to  show  that  | a(2)  +  b|>l  .  Towards  proving  this,  let 
y=sa  +  b»b(a/b  +  l)s=b(x  +  l)  or  x  =  (y/b)-l  and  substituting  this  in  (65) 
yields  the  polynomial 

V(y)  =  y3-3by2  +  (3  +  p)b2y  +  (-  l  -  p  +  q )  b3  . 

Let  Ji,  Y2  and  y3  denote  the  three  real  roots  of  this  cubic  equation.  We  will  show 

that  all  these  roots  have  magnitudes  greater  than  unity.  In  fact,  from  (77)  and  the 
assumptions,  we  have  yi  =  a(1)  +  b>l  and  further  using  (73), 

V(0)  =  (-l-p  +  q)b3>0  .  Since  one  of  its  roots  yi  is  greater  than  unity  and  the 
ordinate-intercept  is  positive,  to  establish  the  above  claim,  it  is  enough  to  show  that 
V(  1  )  >  0  and  V(- 1)  >  0  .  In  that  case,  due  to  its  cubic  nature,  V(y)  must  have  the 
shape  shown  in  Fig.  3(a)  which  corresponds  to  a  strict  Hurwitz  polynomial.  An 
easy  manipulation  shows  that 
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which,  in  particular,  establishes  that  |  y2 1  =  |  a(2)  +  b  |  >  1.  Thus,  p2(z)  =  1  /  (  a(2)  +  bz ) 

also  leads  to  a  b.  r.  function  and  consequently  to  a  solution  to  the  original  problem. 

To  complete  the  proof,  it  remains  to  exhaust  the  remaining  choices  for  b 
and  xj,  viz  ;  (  b  >  1 ,  x1  <  0  ;  b  <  - 1 ,  xx  >  0 ;  b  <  - 1  ,  x1  <  0 ).  Of  these,  for  example, 

when  b>l  and  XjcO  ,  we  have  a(1)<0  and  a(2)>0  and  a(3)>0  and  hence 
|a(2)  +  b|>l  is  automatic.  In  that  case,  we  need  j a(2)  —  b | >  1  and  letting 
z  =  a-  b  =  b(x-l)orx  =  (z/b)  +  l  ,  leads  to  the  new  cubic  equation 

W(z)  =  z3  +  3bz2  +  (3  +  p)b2z  +  (l+p  +  q)b3  . 


Once  again  since  -  a(1)  -  b  <  - 1  and  W(0)  <  0  ,  to  prove  |  Z2  j  >  1  ,  it  is  enough  to 
show  that  W(  1 )  <  0  and  W(- 1)  <  0  .  In  fact, 


W(l)  =  -(b  +  1 ) 


-1 +b  <  0 


and 


W(-l)  =  -( b-1 )( — -1  +  b 


<  0 


and  this  situation  corresponds  to  that  in  Fig.  3(b).  Thus,  |  ^  |  =  |  a(2)  —  b  |  >  1  ,  etcSl8) 

The  remaining  case  can  be  dealt  with  in  a  similar  manner.  Notice  that,  to 
determine  the  bounded  reality  for  P2^z)  ,  we  have  made  use  of  the  previously 
established  fact  that  |  a(1)  ±  b  |  >  1 . 

The  general  case  involving  complex  correlations  can  be  dealt  with  in  a 
similar  manner.  In  that  case,  to  exhibit  an  admissible  solution,  it  is  enough  to 


(18)  Interestingly,  the  above  arguments  also  establish  that  |  a(3)  ±  b|  >  1  .  However,  this  is 
irrelevant,  since  |  x3  |  <  1  . 


demonstrate  the  existence  of  a  bounded  function  (p(z)  is  bounded  if  |  p(z)  |  £1  in 
I  z  |  <  1  ). 

Figures  4-6  together  with  Tables  1-3  show  details  of  maximum  entropy 
extension  and  the  two-step  predictors  discussed  above  in  three  different  situations. 
In  Fig.  4,  the  original  spectrum  corresponds  to  an  autoregressive  model  of  order  6 
and  the  number  of  known  Fourier  coefficients  is  taken  to  be  7.  As  pointed  out  in 
the  Conclusions  (Section  IV),  the  maximum  entropy  extension  results  in  the 
same  autoregressive  form  as  the  original  one  and  the  maximally  robust  two-step 
predictors  generate  stable  ARMA(6,  1)  forms  with  spectra  as  shown  in  Fig.  4  (c). 
The  dotted  curve  represents  the  spectrum  corresponding  to  X!  and  the  solid  curve 
that  due  to  x2.  The  corresponding  filter  coefficients  are  tabulated  in  Table  1.  In 
Figs.  5  and  6,  the  original  spectra  correspond  to  those  of  ARMA(8, 1)  and  ARMA(6, 
3)  processes  respectively.  In  both  cases  the  number  of  known  Fourier  coefficients 
is  taken  to  be  p  +  1,  where  p  corresponds  to  the  respective  denominator  degree. 
The  associated  maximum  entropy  extension  as  well  as  the  best  two-step  predictor 
extensions  are  shown  in  Fig.  5  (b)-(c),  Fig.  6  (b)-(c),  and  the  corresponding  filter 
coefficients  are  tabulated  in  Table  2  and  Table  3,  respectively.  Notice  that  in  all 
these  cases,  the  dotted  curve  appears  to  be  closer’  to  the  maximum  entropy 
extension,  while  the  solid  one  tends  to  'follow'  the  original  one.  This  is  not 
surprising,  since,  from  the  discussions  in  Section  IV,  the  entropy  corresponding 
to  the  dotted  curve  is  greater  than  that  of  the  solid  one  and  consequently  the  dotted 
curve  is  'closer'  to  the  maximum  entropy  extension.  However,  since  these 
extensions  only  possess  a  single  zero,  they  cannot  truly  approximate  spectra 
containing  a  larger  number  of  zeros. 

To  summarize,  through  a  constructive  procedure  we  have  demonstrated 
the  existence  of  two  Wiener  factors  that  maximize  the  minimum  mean  square 
error  associated  with  two-step  predictors  that  are  compatible  with  the  given 
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correlations  r0,  rj, . . .  rn.  They  turn  out  to  be  stable(19)  ARMA  (n,  1)  filters  given 
by 


B^z)  = 


a  +  Pz 


I  gkZk 

k-0 


where 


with 


and 


a  =  ^[V(a  +  b)2-l  +  V(a-b)2-l  ], 
P  »  ^-[V(a  +  b)2  -1  -  V(a  - b)2 - 1  ] , 

A 


a  =  -  2  b  R  cos 

9  ' 

q 

or  -  2  b  R  cos  [  2.  +  — 

U 

l3  3 

7 


(79) 


and  gfe,  k  =  0  -»  n,  as  given  in  (59).  Further,  R  and  <p  are  as  in  (70)  and  (72), 
respectively.  As  remarked  earlier,  the  two  choices  for  the  parameter  a  in  (79)  give 
rise  to  two  bounded-real  functions  and  two  admissible  Wiener  factors. 


IV.  Conclusions 

In  this  paper,  through  a  constructive  procedure,  we  have  demonstrated  the 
existence  of  two  admissible  power  spectra  with  the  property  that  the  associated 
Wiener  factors  maximize  the  minimum  mean  square  error  among  the  class  of  all 
two-step  predictors,  that  are  compatible  with  the  given  (n+1)  correlations  of  a 
stationary  stochastic  process.  Maximization  of  the  two-step  prediction  error  is 


(19)  Note  that  Dn(z)  is  strict  Hurwitz  so  long  as  p(z)  is  b.  r. 
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shown  to  result  in  two  unique  and  stable  ARMA  (n,  1)  filters  with  details  as  given 
in  section  III.  This  is  a  true  generalization  of  Burg's  maximum  entropy 
extension  which  results  in  a  stable  AR  (n)  filter. 

As  remarked  before,  from  a  geometric  viewpoint,  the  maximum  entropy 
extension  is  also  maximally  robust  because  the  new  estimates  so  obtained  are 
always  at  the  centers  of  the  admissible  circles  (see  (25)).  As  a  result,  in  that  case 
the  radius  b0  of  'the  final  circle1  has  the  maximum  possible  value  • 

Naturally,  any  other  extension  cannot  have  all  its  new  estimates  at  the  centers  of 
admissible  circles  and  from  (20)  -  (22),  it  now  follows  that  the  radius  of  'the  final 
circle'  for  any  other  admissible  extension  will  be  strictly  less  than  that 
corresponding  to  the  maximum  entropy  extension.  Thus,  the  radius  of  'the  final 
circle'  may  be  taken  as  a  measure  of  the  robustness  for  the  extension  under 
consideration.  In  the  case  of  ARMA(n,  1)  extension,  from  (21)  and  (60),  the  'final 
radius'  is  given  by 


b0 


(V(a  +  b)2-l  +  V(a  -  b)2  -1 


2a 


(80) 


where  |a|>|b|  +  l  and  |b|>l.  Since  the  radii  of  admissible  circles  form  a 
monotone  nonincreasing  sequence,  this  implies  that  gill  ARMA  (n,  1)  extensions 
have  a  'final  radius',  which  is  at  least  70%  of  the  maximum  possible  value.  Since 
the  specific  ARMA  (n,  1)  extension  discussed  in  Section  III  maximizes  the  two- 
step  prediction  error  ( |  bo|2  + 1  bi  |2  )  ,  the  final  radius  in  that  case  should  possess  a 
tighter  lower  bound. 

Given  a  finite  number  of  covariances,  so  far  we  have  demonstrated  two 
admissible  Wiener  factors  that  maximize  the  two-step  prediction  error  by 
exhibiting  two  b.  r.  functions  from  the  real  roots  of  a  cubic  equation.  The  special 
structure  of  the  numerator  part  of  the  Wiener  factor  in  (48)  together  with  its 


-65- 


general  form  in  (33)  gave  rise  to  the  above  cubic  constraint  that  must  be  satisfied 
by  all  admissible  b.  r.  functions.  The  particular  solutions  presented  in  Section  HI 
involve  two  of  its  roots  -  the  largest  ones  in  magnitude  -  and  since  there  are  two 
admissible  candidates,  this  raises  the  natural  question  :  Does  the  particular 
solution  given  above  using  the  largest  root  have  any  unique  property? 

Interestingly,  in  such  a  situation,  the  solution  corresponding  to  the  largest 
root  possesses  greater  entropy  (  i.e.,  larger  'final  radius’)  and  in  that  sense  it  is 
most  robust  among  the  two  solutions.  This  follows  easily  by  noticing  that  the 
entropy  H(a)  of  an  ARMA  (n,  1)  admissible  extension  is  given  by  H(a)  -  In  b^  ,  with 
b0  as  in  (80)  where  a  and  b  as  a  pair  satisfy  |a|>|l>|  +  X,  |b|»|  a<> / a„ | >  1 .  In  that 
case,  it  is  easy  to  show  that  9H  /  9a  >  0  and  further  a  =  <*>  is  its  only  stationary 
point.  Hence,  H(a)  is  a  monotone  increasing  function  of  a.  Referring  back  to  (76), 
since  a(1)  =  xl  b  is  larger  than  a(2)  =  x2  b  in  magnitude,  we  have  H(  a(1) )  >  H(  a(2) ) 
proving  our  assertion.  In  this  sense,  the  ARMA(n,  1)  extension  characterized  by 
the  largest  root  in  Section  III  is  unique  :  it  maximizes  the  two-step  prediction 
error  and  possesses  maximum  entropy  between  the  two  admissible  solutions. 

The  uniqueness  of  the  above  ARMA  (n,  1)  should  not  be  confused  with  other 
admissible  ARMA  (n,  1)  extensions.  In  fact,  from  (51),  with  b  as  given  by  (54)  and, 
for  example,  letting  |a|>|b|  +  l  results  in  a  bounded  real  function  p(z)  that 
generates  an  admissible  ARMA  (n,  1)  extension.  Thus,  there  is  an  infinite 
number  of  admissible  ARMA  (n,  1)  extensions  that  match  the  first  (n+1) 
correlations,  and  the  ;  articular  one  described  above  is  distinguished  by  the  fact 
that  it  is  most  robust  (maximizes  the  minimum  mean  square  two-step  prediction 
error  and  has  maximum  entropy)  among  all  such  ARMA  (n,  1)  admissible 
extensions. 

The  existence  of  a  single  AR  (n)  type  of  extension  as  well  as  the  abundance 
of  admissible  ARMA  (n,  1)  extensions  that  match  the  given  (n+1)  correlations 
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follow  from  a  general  result  on  Pad6  approximations  [19].  More  generally,  given  n 
=  (p  +  q  +  1)  correlations,  an  ARMA  (p,  q)  admissible  extension,  if  it  exists,  is 
unique,  i.e.,  there  can  be  only  one  admissible  extension  of  the  ARMA  (p,  q)  type 
that  matches  the  first  (p  +  q  +  1)  correlations.  Of  course,  if  n  <  (p  +  q  +  1),  then  the 
admissible  extensions  need  not  be  unique. 

Finally,  given  a  set  of  (n+1)  correlations,  the  question  of  maximizing  the 
entropy  among  all  admissible  ARMA  (n,  1 )  type  extensions  is  also  interesting  by 
itself.  From  (31)  and  (80),  this  is  equivalent  to  maximizing  the  above  b0with 
respect  to  the  variable  a.  As  pointed  out  before,  9b0  /  9a  =  0  results  in  the  solution  a 
=  °o  and  this  does  not  lead  to  a  realizable  b.r.  solution  p(z).  Thus,  given  (n+1) 
autocorrelations,  an  ARMA  (n,  1)  type  admissible  extension  that  maximizes 
entropy  does  not  exist. 

Further,  the  natural  generalization  to  k-step  predictors  (k  >  2)  that 
maximize  the  corresponding  minimum  mean  square  error  is  shown  to  result  in  a 
structured  ARMA  (n,  k-1)  filter  in  Appendix  A.  In  general,  to  identify  an 
ARMA(p,  q)  system,  (p+q+1)  correlations  are  needed.  But  if  fewer  correlations  are 
to  be  used  to  identify  the  same  system,  then  this  approach  leads  to  a  feasible 
solution.  In  that  case,  begining  with  (p+1)  correlations,  maximization  of  (q+1)- 
step  prediction  error  leads  to  an  ARMA(p,  q)  solution.  The  parametrization  of 
such  filters  is,  of  course,  much  more  complicated. 

Appendix.  A 

Maximization  of  the  k-step  Prediction  Error 

Given  a  finite  set  of  (n+1)  autocorrelations  r0,  rj,  r2,  . .  .  rn,  the  problem  here 
i3  to  find  the  extension  that  maximizes  the  minimum  mean  square  error 
associated  with  the  k-step  predictor.  Using  (29)  and  (30),  the  k-step  minimum 
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mean  square  prediction  error  is  given  by 

k-l 


Pk  =  X  lbi!2  =  ed°  ( 1  + 1  cxi  |2  + 1  o&2 12  +  -  •  *  + 1  Ok-i  |2  )  (A.1) 


with  d0  as  defined  in  (26).  Maximization  of  with  respect  to  the  free  parameters 
^n+l'  rn+2>  •  ■  •  f®*^ds  to 

apk  ,  , 

- —  =  0  ,  |  m  |  >  n 


and  using  (A.1),  this  is  equivalent  to 

/k-l  \  ^  /k-l  \ 

(  2  KP  J  +  g^l  S  |“i|2)  =  °-  |m|>n. 

C'rm  \i*0  J  \i  =  0  J 

We  will  show  that  the  left  hand  side  of  the  above  expression  reduces  to 


(A.2) 


k-l 

1 

2 

S(0) 

-n 

i»0 

(A.3) 


and  this  in  turn  implies  that  the  periodic  function  given  by  (A.3)  is  identically 
equal  to  zero  for  all  |  m  |  >  n.  As  a  result,  the  nonnegative  periodic  function 

X  k-l  2 

—  Y  a,  ei®  must  truncate  after  n-terms  and  hence 
S(9)  i-o 


k-l 

^  n 

n 

2  ^ 

=  2  fie*9  = 

X  Si  e*0 

i-0 

i*-n 

i  =  0 

k-l  2 

SC0)  =  - j  =  |  Bk<  eJ0 )  [2 

I  gi e*9 

i-0 


(A.4) 
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where 


k-l 

X  aizl 

Bk(z)  =  ^ -  -  ARMA(n,  k-l) . 

X  SiZ1 

i-0 


(A.5) 


Thus,  to  complete  the  proof,  it  is  enough  to  establish  the  equality  between  the  left- 
hand  side  of  (A.2)  and  (A.3).  Towards  this,  from  (26),  we  obtain 


ddo 

arm 


e->m0  d0 


(A.  6) 


and  by  direct  expansion 


k-l  2  k-l  k-l  p-1 

X  =  X  |ai|2+  X  X  2  Re  (  Op  ejpea*  e_jq0 )  . 

«  A  A  A  A  '  ‘  ' 


(A.  7) 


p*0 qaO 


Substituting  (A.6)  in  (A.2)  and  (A.7)  in  (A.3),  the  desired  equality  condition  between 
(A.2)  and  (A.3)  reduces  to 


k~1  ^ 

X  — 

pTo  arm 


k-l  p-1 


— |«p|2  =  —  2  X  2Kc(apei>*-a,‘e-«9)ejm9de,  |m|>n.(A.8) 

rm  2  7C  p  ■  0  q  =  0 


However, 


p=o  a^m 


«p|2=  X  ^-ke*»|2=  2  2 Jto /«,.*•.  M 

P-0arm  ^  pTo  \  arm ' 


and  hence  with  the  above  expression  in  (A. 8),  a  term  by  term  comparison 
simplifies  (A.  8)  into 
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i-kV**)  =  -M  I  «,VW»  e^dS, 
2x  S(0)U-o  / 


(A.9) 


where  m  >  n  and  0  <  p  <,  k  - 1 


We  will  establish  (A.9)  through  an  induction  argument  on  p.  To  avoid 
cumbersome  notations,  we  will  assume  (A.9)  is  true  up  to  p  =  k  - 1  and  will  extent 
that  result  to  p  =  k.  Clearly  (A.9)  is  well  known  for  k  =  1  and  is  proved  for  k  =  2  in 
Section  III.  To  make  further  progress  for  k  >  2,  notice  that  with  ock  and  dk  as 
defined  in  (29),  straightforward  algebra  gives  the  useful  identity  [1,10] 

ak  ss  dx  ak-1  +  2  d2  ak_2  +  •  •  •  -Kk  - 1)  dk_!  ax  +  k  dk  Oq  ) ,  k  >  1 .  (A.10) 

K 

and  hence 


a<*k_l  y  *  V  i*  aock-i 
3rm  "  k  itl  drm  ^  +  iti  1  drm 


k  r 

1  Y1  *  1  i 

k  iti  2k  £ 

/  -n 


-X 


where  we  have  made  use  of  (26)  and  the  identity  <Xq  =  1 .  By  assumption,  since  (A.9) 
is  true  for  p  =  k  - 1 ,  substituting  that  into  the  above  expression,  we  obtain 


9ak  e 


_1_  A(0) 

2k  S(0) 

J  -n 


ijme  d0 


(A.11) 


where 


k-i-l 


A(0)  =  1+f  £  i  “k-i e  j(k_j)0  +  d*  e*9  £  a*  e'jq 

k!  i-l  q-0 


(A.12) 
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Expanding  the  terms  in  (A.12)  and  rearranging  them  with  the  help  of  (A.10),  it 

k 

simplifies  into  £  “i®-*9.  This  establishes  the  identity  (A.9)  for  p  =  k,  and  completes 
1*0 

the  proof. 

Thus,  given  r0,  ri,  r2,  •  •  •  rn ,  the  Wiener  factor  that  maximizes  the  k-step 
minimum  mean  square  prediction  error,  has  the  form 

k-l 

S  “rn2” 

Bk(z)  =  ^ -  ~  ARMA(n,  k-l) . 

i  gm*m 

m*0 

Notice  that  the  coefficients  in  the  numerator  factor  here  are  precisely  those  given 
by  (29)  and  moreover,  from  (33),  this  extension  should  also  follow  for  a  specific 
choice  of  rational  p(z)  with  numerator  degree  (k-2)  and  denominator  degree  (k-l). 
The  denominator  degree  requirement  for  Bjc(z)  together  with  the  above 
restrictions  on  its  numerator  factor  completely  specifies  the  b.  r.  function  p(z)  and 
hence  the  Wiener  factor.  For  the  existence  of  such  a  Wiener  factor,  the  set  of 
equations  so  obtained  must  yield  a  b.  r.  solution.  Whether  this  is  always  possible 
for  k  >  3,  is  still  an  open  question. 


Appendix  B 


In  this  appendix,  we  will  show  that  the  discriminant  D  in  (69)  is  always 
negative.  To  prove  this,  with  (67)  -  (68)  in  (69)  and  letting  aj/a^  =  y  ,  the 
discriminant  D  simplifies  into 


q 

|2 

4- 

P 

2 

3  ~  1 

3 
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’i(1+yy'Hy'^[y3+6(1+p)y2+8(1+p)T  y>0-  ®a) 

l  y  y 

Since  |b|>l,  we  have  <1,  which  gives  ^4  <^2  with  this  in  the  above 


expression,  it  simplifies  into 


D  s  ly-iPy-^[y3+6(1+^|y2+8(1+^)1  ■  -i^y-^f  ■  «® 


where,  by  definition 


A(y)  =  y3  +  6  (l  +  y2  - 15  y  +  8  (l  +  . 


(B.3) 


Clearly,  to  establish  D  <  0,  it  is  enough  to  show  that  A(y)  >  0  for  all  y  >  0  and  |  b  |  >  1 . 


Let  (3  =  2  |l  +  ,  then  I  b  I  >  1  implies  2  <  (3  <  4  and  (B.3)  simplifies  into 

A(y)  m  y3  +  3  (3  y2  - 15  y  +  (33  . 


(B.4) 


The  desired  result  follows  if  Amin  >  0  for  y  >  0  and  2  <  |3  <  4.  To  find  its  minimum, 


note  that 


=  3y2  +  6  Py-15  =  0 

dy 


yields  the  positive  solution 


y0  =  -p  +  V|32  +  5 


and  the  minimum  value  for  (B.4)  is  given  by 


A(y0)  =  3(33  -  2  ((32  +  5)  V  p2  +  5  +  15(3  =  Ao({3)  . 


(B.5) 


As  Fig.  7  shows  A(y0)  >  0  for  2  <  J3  <  4  and  this  completes  the  proof. 
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Appendix  C 
Uniqueness  of  p(z) 


For  uniqueness  of  p(z)  in  (51),  it  remains  to  show  that  other  feasible  forms 
such  as 


p(z)  = 


1  +  djz  +  d2z2  +  •  •  •  +  duj.iZ131-1 
a  +  bz  +  cxz2  +  c2z3  +  •  •  •  +  Cnj^z™ 


m  >  1 


(C.1) 


cannot  give  rise  to  bounded  real  functions  of  degree  one  in  (48),  unless  q  =  d*  =  0,  i 
=  1,2,...  m-1.  To  prove  this,  consider  first 

p(z)  =  i  +  &  .  (C.2) 

a  +  bz  +  cz* 


From  (33),  since  B2(z)  =  r(z)/Dn(z),  the  degree  one  requirement  in  the  numerator  of 
the  Wiener  factor  B2(z)  in  (48),  translates  into  the  same  requirement  for  the 
numerator  of  Hz).  Thus  with  (C.2)  in  (34),  we  obtain 

(a  +  bz  +  cz2)  (a  +  bz  +  cz2)*  -  (1  +  dz)  (1  +  dz)*  =  (a  +  (3z)  (a  +  (3z)* 
which  is  the  same  as 


ac  (z2  +  z~2)  +  (ab  +  be  -  d)  (z  +  z  x)  +  a2  +  b2  +  c2  -  d2  - 1  =  ap  (z  +  z_1)  +  a2  +  p2 . 


Comparing  the  equal  powers  of  z  on  both  sides,  we  must  have  either  a  =  0  or  c  =  0. 
If  a  =  0,  then 


p(z)  = 


1  +  dz 
z  (b  +  cz) 


and  it  is  not  a  b.  r.  function  because  of  the  pole  at  the  origin.  On  the  other  hand,  if 
c  =  0,  then 


p(z)  = 


1  -»•  dz 
a  +  bz 
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and  substituting  this  in  (35),  we  obtain 


-  da<)Zn+2  +  (ban  ~  &o  +  da1)zn+1  + 
a  +  bz 


Dn(z)  = 


Since  Dn(z)  is  of  degree  n,  necessarily  dag  =  0  =*■  d  =  0,  since  ao  *  0  and  this  leads 
to 

1 


p(z)  = 


a  +  bz 


(C.3) 


Arguing  in  a  similar  manner,  starting  with  any  m  >  2,  (C.l)  can  be  shown  to 
reduce  to  (C.3).  This  proves  our  claim. 
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Fig.  2  Structure  of  the  cubic  function  x3  +  px  +  q,  p  <  0  (see  (65)-(68)). 
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(b)  One  Step  Predictor(Maximum  Entropy  Estimator)  AR(6)  -  B^z) 


Fig.  4  One  and  two  step  predictors  that  maximize  the  minimum  mean  square  error. 
Original  spectrum  corresponds  to  an  Autoregressive  model  of  order  6  given  by 


1.340  -  3.044z  +  4.599z2  -  4.825z3  +  4.171z4  -  2J04z5  +  z 
The  number  of  Fourier  coefficients  known  is  taken  to  be  7. 
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Table  1 


One-Step  and  Two-Step  Predictors 
that  Maximize  the  Minimum  Mean  Square  Error 

Q(z)  Ql(z)  Q2^Z) 

B(z)  ■  rtf ;  B,(z)  -  w ;  b2(2)  -  —  . 

B(z) ,  Bt(z)  and  B2(z)  represent  the  original,  the  one-step 
and  two-step  Wiener  factors,  respectively. 

n  =  No.  of  Fourier  coefficients  used  for  (b)  and  (c). 


Fig.No. 

n 

Model 

Wiener  Factor  Coefficients 

a 

(original) 

AR(6) 

Q(z)  =  1 

P(z)  =  1.340  -  3.044z  +  4.599z2  -  4.825z3  + 
4.l71z4  -  2.504z5  +  z 

Fig.4 

b 

7 

AR(6) 

Oft)  =  i 

Pt(z)  =  1.340  -  3.044Z  +  4.599z2  -  4.825z3  + 

4.171z4  -  2.504z5  +  z 

c 

(dotted) 

ARMA(6,1) 

Q2(z)  =  -  0.895  +  0.322z 

P2(z)  =  1.243  -  3.009z  +  4.632z2  -  4.948z3  + 
4.25 0z4  -  2.553z5  +  z 

c 

(solid) 

ARMA(6,1) 

Q;(z)  =  1.044  +  0.4 16z 

P;(z)  =  1.464  -  3.089z  +  4.557z2  -  4.680z3  + 
4.070z4  -  2.442z5  +  z 
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(b)  One  Step  Predictor( Maximum  Entropy  Estimator)  AR(8)  -  Bt(z) 


Fig.  5  One  and  two  step  predictors  that  maximize  the  minimum  mean  square  error. 
Original  spectrum  corresponds  to  an  ARMA(8,1)  given  by 


3(z)  = 


_ 1.1  +  z _ 

1.267  -  0.635z  +  r  3 66z  -  0.664z3  +  1.425z4  -  0.626z5  +  0.770z6  -  0.53  lz7  +  z8 


The  number  of  Fourier  coefficients  known  is  taken  to  be  9. 
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Table  2 


One-Step  and  Two-Step  Predictors 
that  Maximize  the  Minimum  Mean  Square  Error 

Q(z)  Ql(Z)  Q2(z) 

BW  ■  m  '•  B*(z) - m ;  Bi(z) = w  • 

B(z),  B^z)  and  B2(z)  represent  the  original,  the  one-step 
and  two-step  Wiener  factors,  respectively. 

n  =  No.  of  Fourier  coefficients  used  for  (b)  and  (c). 


Model 

Wiener  Factor  Coefficients 

ARMA(8,1) 

Q(z)  =  1.1  +  z 

P(z)  =  1.267  -  0.635z  +  0.866z2  -  0.664z3  + 
1.425z4  -  0.626z5  +  0.770z6  -  0.53  lz7  +  zS 

AR(8) 

Qi(z)  =  i 

Px(z)  =  1.411  -  1.1 15z  +  1.367z2  -  1.132z3  + 
1.830z4  -  1.007z5  +  0.976z6  -  0.663z7  +  z8 

ARMA(8,1) 

Q2(z)  =  -  1.004  +  0.6  lOz 

P2(z)  =  1.563  -  1.623z  +  1.898z2  -  1.627z3  + 
2.258z4  -  1.4 10z5  +  1.1 94z6  -  0.663z7  +  z8 

ARMA(8,1) 

Q2(z)  =  0.766  +  0.594z 

P2(z)  =  1.272  -  0.65  lz  +  0.883z2  -  0.680z3  + 
1.439z4  -  0.639z5  +  0.777z6  -  0.533z7  +  z8 

(a)  Original  Spectrum  ARMA(6,3)  ~  B(z) 


_ 

\ 

\ 

\. 

(b)  One  Step  Predictor(Maximum  Entropy  Estimator)  AR(6)  ~  B.(z) 


(c)  Two  Step  Predictors  ARM A( 6,1)  -  B,(z) 

Fig.  6  One  and  two  step  predictors  that  maximize  the  minimum  mean  square  error. 
Original  spectrum  corresponds  to  an  ARMA(6,3)  given  by 


-  1317  -  1.089z  +  1.181z2  +  z 


1.340  -  1.201z  +  0.96 lz2  -  1.208z3  +  0.871z4  -  0.988z5  +  z6 


The  number  of  Fourier  coefficients  known  is  taken  to  be  7. 
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Table  3 


One-Step  and  Two-Step  Predictors 
that  Maximize  the  Minimum  Mean  Square  Error 


B(z)  =  ;  B.(z)  = 

V  ;  P(z)  ’  lK  '  Px(z) 


Qi(z>  Q2(z) 

;  b2(z)  = 


P2(z) 


B(z) ,  Bt(z)  and  B2(z)  represent  the  original,  the  one-step 
arid  two-step  Wiener  factors,  respectively. 


n  =  No.  of  Fourier  coefficients  used  for  (b)  and  (c). 


Fig.No. 

n 

Model 

Wiener  Factor  Coefficients 

Fig.6 

a 

(original) 

7 

ARMA(6,3) 

Q(z)  =  -  1.317  -  1.089z  +  1.181z2  +  z 

P(z)  =  1.340  -  1.201z  +  0.96 lz2  -  1.208z3  + 
0.87 lz4  -  0.988z5  +  z 

b 

AR(6) 

Qi(z)  =  i 

Pr(z)  =  1.432  -  1.462z  +  1.645z2  -  1.660z3  + 
1.279z4  -  1.129z5  +  z 

c 

(dotted) 

ARMA(6,1) 

Q2(z)  =  -  0.879  +  0.485 z 

P2(z)  =  1.363  -  1.720z  +  1.867z2  -  1.917z3  + 
1.442z4  -  1.133  +  z 

c 

(solid) 

ARM  A(  6,1) 

Q2(z)  =  0.934  +  0.649z 

P2(z)  =  1.520  -  1.133z  +  1.36  lz2  -  1.333z3  + 
1.072z4  -  1.124z5  +  z6 
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